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INTRODUCTION 


Since the recent appearance of the branch of structural optimization which 
is concerned with constraints of a dynamic (natural frequencies of vibration) or 
aeroelastic nature*, two trends have developed, which represent two complemen- 
tary methods of solution. 

The first of those makes use of mathematical programming techniques. 
The structure to be optimized is either a structure with discrete characteristics, 
such as a stiffened beam, or a continuous one which is first discretized: this is 

the case of a beam broken into finite elements. This approach, the older, is the 

2 3 4 5 6 

one taken by Haug, et al. ’ , Fox and Kapoor , Rubin , and Turner , to name 

only a few. Its advantages are obvious: complex structures can be approximated 
this way with good accuracy, then subject to optimization techniques. The weak- 
ness of the method, however, lies in the fact that the optimization is done after 
the discretization has been made, so that it is impossible to know how close the 
discrete optimum found is from the true one, and even in some extreme cases it 
makes no sense at all. ^ In other words, there is no way of getting the optimum 
within a given accuracy. 


*A small paragraph was devoted to the subject under the heading "dynamic com- 
pliance" in the excellent comprehensive review of the field of Structural Optimi- 
zation done by Sheu and Prager^ less than three years ago. The number of 
articles belonging to this classification which have appeared since is impressive, 
and the subject alone now merits far more extensive review. 

^The only way to decide where to stop iterating in the optimization process is when 
it appears to converge; however, it has been frequently proven that, depending on 
the method chosen to arrive at the optimum, the process might pseudo-converge 
to a solution which is totally irrelevant and vanishes when another optimization 
method is tried. 


Another approach, more recent, is that of the classical calculus of 

variations whose foundations date back to the time of Euler and Lagrange and has 

raised considerable interest ever since: it is applicable to simple, continuous 

7 

structures. It has been widely used by Turner for his now classical bar, by 

8 9 

Keller and Niordson for the problem of the tallest column, by Niordson for 

the optimization of a beam for a fixed first frequency of vibration, by Olhoff 10 

for the optimal design of vibrating circular plates. Variational principles were 

11 12 

also used by Prager and Taylor and Taylor and applied to continuous 
systems. 


In the spirit of this approach an extension was made by the use of the 

techniques of optimal control theory, which are themselves ramifications of 

the classical calculus of variations; its principles are expressed in a different 

form, more suitable to the kind of problems encountered in the theory of control. 

For the structural applications, the constraints are put under a form analogous 

to that of control constraints , thus imposing the distinction between state and 

control variables. The necessary conditions are then derived in the form of 

a system of differential equations, yielding a two-point boundary -value problem. 

Exact solutions have been found in a few rewarding cases, and numerical methods 

derived in control theory, requiring discretization but only at this step, have 

been applied. Accuracy can then be fixed in advance, the limitations of this 

method now being due to the complexity of the structure. In other words, this is 

the opposite approach from the one discussed first. This method was first 

13 

introducted by Ashley and McIntosh for the optimization of structures subject 

14 

to aeroelastic constraints, and applied by Armand and Vitte to a few simple 

15 

cases and by Weisshaar to the panel flutter problem started by Turner. 

A recent and detailed account of the development of the field of structural 
— especially dynamic or aeroelastic — optimization will be found in some of the 
references above, especially ^ and in references"*^’ 

The structures optimized up to now were "one-dimensional"; this means 
that, due to their nature, the problems to be investigated led to constraints 
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expressed in the form of ordinary differential equations in one independent 
spatial variable. A beam problem is obviously one-dimensional; so are optimi- 
zation problems dealing with circular plates when the imposed static or dynamic 
constraints have the property of axisymmetry. The problems of dynamic structural 

optimization where the constraints are in the form of partial differential equations 

have hardly been investigated up-to-date. The only tentative efforts known to the 

18 19 

author are those of Johnson and Haug who investigated separately the problem 
of the optimization of a plate for a fixed first frequency of vibration. Both of them 
make use of the method of discretization described at the beginning of this chapter, 
and this is why, unfortunately, it is impossible to give an estimate of the pre- 
cision which guarantees the validity of their optimal solution. Moreover, 

Johnson' s square plate is divided into only 25 elements, the symmetry of the 
problem further reducing the accuracy to only 6 elements. Haug, on the other 
hand, presents an adaptation of the powerful method of steepest descent developed 
in optimal control theory to two-dimensional problems , which seems very 
promising and should lead to the numerical solution of number of such problems. 

The present work tries to present an adaptation to two-dimensional 
structures of the methods of optimal control theory already applied to one- 
dimensional structural optimization, in the spirit of the previous work done in 
this area^’ ^ The corresponding branch of Optimal Control Theory is 

known as optimal control of distributed-parameter systems and is, we might 
say, a brand new field of investigation. The first ever made, proposing a 

definition and an evaluation of the task at hand, is contained in a Russian paper 
20 

dated 1960 The systems to be considered are characterized by either partial 
differential equations or integral equations to be satisfied on DxT, where D is 
some spatial domain contained in the Euclidean space E n for n ^ 1 and T is 
a time domain. In the case of differential equations, there are also boundary 
conditions to be satisfied. The behavior of the system can be influenced either 
by mechanisms acting throughout D or at the boundary of D, or both. There is 
given some functional dependent on the state of the system or its boundary 
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conditions or on its control, or some combination of these, and two types of 
questions may be asked: 

i) What should the available controls be to minimize the given functional 
(open- loop control) ? 

ii) What should the functional relationship between the control and the 
system state be in order that the functional be minimized (closed-loop control) ? 

Needless to say, optimal control of distributed-parameter systems is 

much more difficult than for lumped-parameter systems, and it is still a growing 

field of investigation. An excellent survey containing an exhaustive list of 

21 

references is due to Robinson . In his paper, both a historical and technical 

presentation of the research in this domain is done. Two textbooks have appeared 

22 

recently, one by the pioneer Butkovskiy presenting a unifying treatment of most 

of the literature that appeared previously in Russian journals and of great interest 

23 

although perhaps too broad in scope; the other one is by Lions and is mostly 
concerned with the mathematical aspect of some particular problems. 

The overwhelming majority of research up-to-date has been devoted to 
linear problems. Although our scope will be more limited than the general goal 
of optimal control of distributed-parameter systems, we shall never encounter 
linear problems, and we will have to derive a more general approach applicable 
to non-linear constraints. 

We first have to define the kind of problems that we will encounter for our 
purpose, which remains structural optimization: the domain DxT will simply be 
the domain occupied in space by the two-dimensional structure which we consider, 
referred to a set of orthonormal axes Ox, Oy. The constraint will take the form 
of a partial differential equation to be satisfied inside the domain, with adequate 
boundary conditions provided by the physical properties of the structure. For the 
common case of a minimum -mass problem, the functional to be minimized will be 
the surface integral of the thickness over the domain covered by the structure. 
Generally, the deflection (in-plane or out-of-plane) will play the role of a state 
variable, whereas the thickness will be in some cases the control variable. 
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A problem of this type is a special case of the more general one stated above, 

namely a Mayer-Bolza problem for multiple integrals. The most authoritative 

treatment of the general Mayer-Bolza problem in two dimensions is that of 
24 

Lur' e , who was the first to suspect the applications this theory might have for 

structural optimization. This paper was later reprinted in a textbook on topics 

25 

in optimization edited by Leitmann with some additions, in particular an 

application to magnetohydrodynamics leading to a control of the bang-bang type, 

itself a reprint of an article published in the Journal of Applied Mathematics and 
26 

Mechanics Lur' e derives the necessary conditions for an optimum for a very 
general Mayer-Bolza type of problem , extending the classical methods of 
calculus of variations to two dimensions. 

27 

In part A of this study, we will extend the methods used by Bryson and Ho 
in their excellent textbook on optimal control theory to the case of two independent 
variables. We will be led to the definition of a scalar two-dimensional Hamiltonian, 
already suggested by Lur' e, and to a set of necessary conditions for an extremal 
in terms of it analogous to the ones he derived by other means. The case of finite 
equality or inequality constraints on the control variables will also be considered. 

Part B presents an application of the theory to the case of the optimization 
of a shear plate: we want to find the thickness distribution such as to minimize 
the mass of this structure, the first frequency of vibration of the homogeneous 
reference plate with constant thickness being required to keep a constant value. 

The shear plate was chosen so as to provide a simple, though non-trivial, partial 
differential equation as a constraint; an analytical solution will be found in the case 
of simply-supported edges. 

Our experience will then be applied to the minimum-mass design of more 
practical structures, under the same constraint as above. The sandwich plate 
will be treated in part C, while part D is devoted to the classical elastic plate. 
Following the steps of the simple shear -plate optimization process, both of these 
optimization problems will reduce into finding the solution of a system of two 
simultaneous partial differential equations for the optimal thickness distribution 
and the displacement mode. The necessary conditions are then in a form showing 
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that they are nothing but the expression of a very general sufficient optimality 
condition stated by Prager for a broad class of structural optimization problems 
to which the considered cases belong; they are therefore necessary and sufficient 
conditions for an extremal. A numerical treatment of this system of partial 
differential equations is presented in Part C for a simply-supported sandwich 
plate. The uniformity of the energy distribution throughout the structure is shown 
in Part E to be a sufficient condition for an extremal of the mass. 


6 



PART A 


N ECESSARY CONDITIONS FOR THE STATIONARY VALUE OF A FUNCTIONAL 
UNDER CONSTRAINTS EXPRESSED AS PARTIAL DIFFERENTIAL EQUATIONS 

1. Statement of the Problem 

Let D denote a closed domain in the plane with a piecewise continuous 
boundary 9D. If the plane is referred to a set of rectangular coordinates (x,y), 
let the boundary curve be represented in parametric form by the equations: 

x = a (ct) 

y = P(o-) 

where a and (3 are piecewise continuous functions of the parameter u possess- 
ing a first derivative in the intervals where they are continuous. 

In this domain, we consider a system of partial differential equations of 
the form: 

9z. 

^-Xg.u^.y)- 0 

9z. 

— Y(z,u;x,y) = 0 (1.1) 

dy i ~ ~ 

i = 1, 2, ... n. 

The X. and Y. functions have to be such that they satisfy the compati- 
bility conditions: 

DX. DY. 

= 0 

Dy Dx 

i = 1,2, . . .n, where the derivatives above are total derivatives taken with respect 
to all the arguments included, i. e. : 

DX. 9X. 9z 9X. 9u 9X. 

l _ i l 

Dy 9z 9y 9u 9y 9y 
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DY. 3Y. 3z 9Y. 9u 9Y. 
i _ i i ^ i 

Dx 9z 9x 9u 3x 3x 


For a system governed by a set of equations in the form (1.1), the vector 
function z = (z , z , . . . , z ) of the arguments x,y describes the mechanical 

\ 2i n 

system itself, and therefore the z. play the role of state variables, whereas 
the functions u = (u ,u , . . . ,u ) of the same arguments represent the distributed 

^ I Z P 

controls . 

Equations (1. 1) represent a standard form of any system of partial 

25 

differential equations: Lur'e calls it the special case of the Pfaffian system. 

Any partial differential equation of order superior to one or any system of such 

equations may be written in the form (1. 1), with the number of dependent variables 

increased as necessary. The technique of decomposition is simply an extension of 

the one used in the case of the transformation of a high order ordinary differential 

equation into a system of first order ones by introducing auxiliary dependent 
14 27 

variables. ’ The only difference resides in the fact that a function in two 
independent variables possesses k+1 derivatives of order k, whereas a func- 
tion in one independent variable possesses only one: therefore the decomposition 
will require the introduction of a much considerable number of dependent variables. 
The compatibility conditions will help reduce this number. 

The procedure is best demonstrated on a simple example, treated below. 

We will then show how the decomposition is done for the general case of a single 
partial differential equation in two dependent variables which we will always 
encounter as a constraint in problems of minimum-mass design of structures 
satisfying conditions on their frequency of vibration or on some aeroelastic 
eigenvalue. Generally, one of the dependent variable will represent the displace- 
ment at a point of the structure, the other representing the thickness distribution. 

For the Helmholtz equation, example also treated by Lur' e in the already 
referenced paper: 


ijL 

2 + 2 + UZ 
3x 9y 


0 
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we find the equivalent system: 


!fi_ 

9x Z 2 


9y Z 3 


9x U 2 


8y u 3 


!!3 

9x U 3 


9y 


"Wl 


Here we had to introduce the new dependent variables z , z ,u ,u , the 

u Z u 

original z being denoted z^ and the original u.u^. The state variables are now 
z 1 ,z_,z , while the controls are u u ,u . 

In the light of the example above, we are now able to explain more precisely 
how the decomposition (1. 1) can be achieved in the case of main concern to us, as 
we will see later, that is of a constraint expressed in the form of one single partial 
differential equation in the two independent variables x and y, and two dependent 
variables w and h of order n ^ 1 and n s 0 respectively in those variables. 

X /L 

In other words, there is present at least one partial derivative of order n^ in the 
first variable called w, and at least one partial derivative of order n 2 in the 
other variable, called h. The most general form for this type of equation will be: 

rt _ XI _ 


f(w, 


9w 9w 9 w 
8 x'3y' 8x 2 


n l 

8 w 9 " L w 


a 1 
o w 


n-i ’ n-,-1 

9x x 9x 1 9y 


9x9y 


V 1 


a n l 

d w 

n i 

3y 


(continued) 
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2 U 2 n 2 

9h 9h 9 h 9 h 9 h 

’ ax ’ 8y ’ 8x 2 .V. W"". . V 1 ’ “2 


n n 

9 2 h 9 2 h 


; x,y) = 0 

9x 9x 9y 9x9y “ 9y (1. 1. si) 


and we will assume that this equation can be solved for one of the highest derivatives 
of w appearing, for example, if 


9 w 


9y 

is present: 


n-i 2 n l 

9 w 9w 9w 9 w 9 w 

= ^ V 9X 2 „ ”1 


9y 


9x 


a n i 
9 w 


n -,- 1 


n l 

9 w 


J • * • J 


9x 9y 9x9y 


■v 1 ’ 


2 n o n 2 

3h % 9_h 9 h 9 h 

' 3X ' 8y ' Sx 2 ’ " " ’ . B 2 ’ “2 _1 

9x 9x 9y 


a\ 


°2 

9 2 h 


n 2 - ^ n 2 
9x9y 9y 


;x,y) (1.1. b) 


n l( n i +1 ) 

We will introduce 1+2+ .. . +n, auxiliary variables z. as 

12 1 


follows: 


w = z. 


!fi_ 

9x Z 2 


(=~) 
1 9x ; 


9y Z 3 


(=~ ) 
' 9y ; 


!!2 

9x Z 4 


9x 2 


9y Z 5 


K dxdy 


(continued) 
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!f3 

9x Z 5 


(-i) 

v 9x9y' 


9y z 6 


2 

«-*f> 

9y 


9z 


(n 1 -l)(n 1 -2) 


+ 1 


9x 


= z 




+ 1 


V 1 

9 w 

(-— — f> 

„ V 1 

9x 


9z 


(n 1 -l)(n 1 -2) 


+1 


= z 


9y 


' n 1 ( n 1 “i) 


+ 2 


V 1 

<-■§ 3E-, 

v n -2 ' 

9x 9y 


9z 


(n 1 -l)(n 1 -2) 


+ 2 


9x 


= z 


n l( n i -1 ) 


+ 2 


a ” 1 " 1 

(-- — ) 

v n -2 ’ 

9x 1 9y 


9z 


(n 1 -l)(n 1 -2) 


+ 2 


9y 


■ = z 


’n 1 ( n 1 -l) 


+ 3 


V 1 

(-1 *— ) 

V 3 2 

9x 9y 


9z 


n 1 (n 1 -l) 


9x 


= z 


m ( v i) 


- 1 


v 1 

9 \v 

(= r — ) 


9x9y 


n i-2 


(continued) 


11 



9z 


n i (n i“ 1 ) 


8 y 


■ = z 


n i (n i + X) 


-V 1 


9y 


V 1 


The (n + 1) u variables are now introduced as follows: 


9z 


n l^ n l -1) 


+ 1 


9x 


= u. 


n i 

(= <LW) 

V 

9x 


9z 




+ 1 


9y 


= u 


n l 

(= — 8 - y > 

k n -1 ’ 

9x 9y 


9z 


9z 


n i (D l +1) 


9x 


n l (n i +1) 


= u 


9y 


= u 


n ^+1 


n i 

C-^r» 


9x9y 


n i 

n l 

9y 


V 1 


We now do the same for the dependent variable h, introducing 


n 2 (n 2 +1) 


auxiliary state variables h^ and n^+1 auxiliary control variables v^ , follow- 
ing an identical procedure. 

Now the given equation will be replaced by the two systems formed above, 
together with the condition: 


V+l cp(z l’ Z 2 ,Z 3” ' ' ’ V(n +1) ’ “l’"' ,u n. , ,h l ,h 2’ V" ’ ’ 


l v 1 


(continued) 
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In the case of a constraint of this type, the role of the control variables is 
played by the u and v variables, while the z and h define the state of the 
mechanical system. The z and h are the state variables z of (1. 1) whereas 
the u and v are the distributed controls u. 

The decomposition of a partial differential equation or of a system of the 
same into a system of the form (1. 1) is of course by no means unique. A question 
of definition of the controls arises too, as can be seen in the following example. 

Consider the partial differential equation, in the two dependent variables 
w and h, that we will encounter as a constraint in part B: 


8 /U \ 9 ,, \ , 2 

5? h 5T ) + 57 (h SF )+khw ‘° 


By the process described above, this equation is equivalent to the system: 


l-H 

N 

II 

■s 


9z l 

dx z 2 

(= — ) 
' 9x' 

9z l 

' 9y 

ay z 3 

9Z 2 

9x U 1 

a 2 

dx 

9Z 2 

ay u 2 

,2 
a w 

8x8y 

9Z 3 

9x U 2 

(- ^ ) 
1 9x8y 

9z 3 

9y u 3 

9y 


(continued) 



ah 

0X 


V 


1 


ah. = 

Sy t 2 


. 2 V 1 Ws 

-a " -v k v — 5 — 


where the role of the control u is actually played by u ,u , v ,v ,u bringing 

f*- 1 1 o 1 z z 

no contribution to the system itself as can be seen from the constraint. Because 
of its special nature, the given equation is also equivalent to the somewhat simpler 
system: 


w = 

N = *2 

3x h 
9y h 

!f2 

ax u i 

3y U 2 
ax u 3 


9Z S .2 

aT = “V k hz i 


where the role of the control is now being played by h,u , u , u , the definition 

1 Z o 

of the u being different than from above. 

We will show in Appendix 1 that, although the two systems and the two 
controls are not the same, the same final result is obtained for the same optimal 
problem they both describe. 
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Let us now go back to the statement of the optimal problem. We now 
have to describe the boundary conditions. 

We assume that the first m s n functions z. are prescribed along 
portions 2 of 9D: 

z.l = z.(a) i=l,2,...m (1.2) 

M2 1 

The simplest Mayer-Bolza problem is now formulated as follows: in a 
suitable class of functions, to which we will come back later, determine the state 
variables z and controls u so as to minimize the functional: 


/ ^ a > da + //'< z,u;x, y)dxdy 


9D 


D 


(1.3) 


subject to side conditions (1. 1) and (1.2). 

The problem may be complicated by considering boundary controls on 
movable boundaries and/or constraints on the state or control variables. The 
case of boundary controls and that of constraints on the state variables will not 
be investigated in the present study; the former can be found, treated in some 
details, in Refs. 24 and 25. However, we will investigate the case of constraints 
(finite equalities or inequalities) on the control variables in section 3 of this 
chapter. 


2. The Necessary Conditions 

Following the methods of optimal control theory in one single variable, we 
adjoin the system of partial differential equations (1. 1) to J with the help of 
vector multiplier functions X(x,y) and y(x,y). 


J = J X(z;u)dc7 + JJ <L(z,u;x, 
9D D ^ 


T 9 ~ 

y) + X (x,y)[X(z,u;x,y) - — ] + 
~ ^ ^ ~ dx 


(continued) 
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3z 


'} 


+ p, (x,y)[Y(z,u;x,y) -— ] } dxdy (2.1) 

r*> r* /v ^ Oy 1 

We define a scalar function, the Hamiltonian, as 

T T 

H(z,u;x,y) = L(z,u;£,y) + A (x,y)X(z,u;x,y) + ll (x,y)Y(z,u;x,y) 


( 2 . 2 ) 


J is rewritten then as: 


3 ^ J' i(zr,cr)da + ff | h ( z , u ; x , 
9D D ^ 

If we add and subtract the quantity: 


T 92 T 9Z> » 

(2.3) 


T T, 

9* 9 Ji 
~ + I z 


y 9x 9y 
to the second integrand, J takes the form: 


T T 

9A 9, 


f Uw)da+jy’Lz,u J x,y) + ^ + f-'jJ\dxdy 
3D D ^ J 


-ff\ 


r\ rri rj rp j 

5 + s)} 


D 


Now, by the use of Green' s formula (*), the third integral becomes a line 
integral as follows: 


M 

D 


3 T 9 , T 1, , 

S ) + ^ ( !i g" fs 


/ T T 
(ji dx - X dy)z 


9D 


(2.3) 


Then the integral (2.3) takes the form: 


(*) which is the counterpart in two dimensions of the classical integration by 
parts in one dimension. 
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J ” [pT'a' (c) - X T P - (a)]zdcr 


3D 

and J becomes: 


J = 


f h(z;a’)-l^ r (o') - g'a' (cr)]:z > do 

3D ^ 


JJ |h(z,u;x, 


y) + 


T T 
3X A 8^ - 

rsJ rw 

3x + 3y 


z > dxdy 


(2.4) 


(2.5) 


Now consider the variation in J due to variations in the control vector 

u(x,y): 


6 J 


f [at 
J _9z 


-X^p' (a) - (i^a' (a) 6 z d a 


3D 



(2.6) 


We now choose the multiplier functions X(x,y) and n(x,y) to cause the 
coefficients of 6z in the surface integral above to vanish, therefore escaping to 
the task of computing the 6z in terms of the variations 6u of the controls: 



3x + 3y 


3H 

3z 


(2.7) 


Along the portions of 3D where the z. are prescribed, 

6 z = 0, i = 1, 2, . . . m. 

Along the other portions and for the z^ , j = m +1 n which are not at 

all prescribed, we ask the relations 
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(2.7a) 


\V (a) - ,/a' (a) = - 

to hold along 3D: they serve as boundary conditions for the system (2. 7). 
(2. 6) then becomes: 


6J = 



dxdy 


For an extremum, 6 J must be zero for arbitrary 6u(x,y); this can only 
happen if: 





( 2 . 8 ) 


Equations (2. 8) are analogous to the control equations derived in one- 
dimensional optimal control theory. 

Equations (2. 7), (2. 7a) and (2. 8) are the Euler-Lagrange equations of the 
classical calculus of variations, for two independent variables. They form a set 
of necessary conditions for an optimum. 

In summary, to find a control vector function u(x,y) that produces a 
stationary value of the performance index: 


/ 

3D 


i(z;cr)dcr + 


// 


L(z, u;x,y)dxdy 


D 


(1.3) 


we must solve the following system of partial differential equations in D: 


3z 


~X(z,u;x,y) 


9z 



(1.1) 

dX an T 

— + — = -(^) 
3x 3y '3z ' 

(2.7) 
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( 2 . 8 ) 


9H 

9u 


= 0 


where the Hamiltonian H is defined as: 


T T 

H(z,u;x,y) =.L(z,u;x,y) + X (x,y)X+ g, (x,y)Y 


r v 


( 2 . 2 ) 


The boundary conditions are: 


z 4 =z * <tr) 


i = 1,2, . . . ,m 


( 1 . 2 ) 


(the 2 are the portions of 9D where the first m variables z^ are prescribed) 


y w - <°) = 


dZ 

9z. ’ 

l 


i = 1,2, ... ,m 


(2. 7a) 


along 9D-S, and: 


dZ 

(O’) - u. a' (o') = - — , j = m+1, . , n 

3 3 j 


(2. 7a) 


along 9D. 

We have n+n+n+p = 3n+p equations for the 3n+p unknowns z,X,^, and 

u. 


3. Eq uality and Inequality Constraints on the Control Variables 

For most of the structural optimization problems that we will consider, 
the role of the control variable will be played by the thickness distribution which 
we wish to render optimal. For obvious physical reasons, we do not want, in 
most cases, for this thickness to be zero along a line interior to the structure 
(leading to the formation of a hinge, subsequently causing the collapse of the 
structure). Such an inconvenience might be taken care of by the adjunction to 
equations (1. 1) and (1.2) of a minimum thickness constraint, which will appear 
as an inequality constraint on the control variables. 
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More generally, suppose we have restrictions imposed on the control 
functions in the form of r finite equalities: 

C(u;x,y) = 0, k=l,2,...r (3.1) 

K ~ 

and s finite inequalities: 

D^(u;x,y) 0, JL = l,2,...,s (3.2) 

We assume that there exists a solution to equations (1. 1), (1. 2) under 
the restrictions (3.1), (3. 2) imposed on the control tunctions. 

A way to transform the inequality constraints (3. 2) into equality con- 
straints is to introduce supplementary artificial controls u* = (u* u* . . . ,u*), 

~ 12 s 

following a traditional technique, by virtue of the following equations: 

D * = D ^0i‘> x - y ) “ = °> X = 1.2 s (3.3) 

and to add these controls to the previous set of u, therefore considering only 

r + s equality constraints with an augmented number of controls. 

24 

This point of view is the one adopted in Lur' e' s paper . We will 
present here another method, again directly inspired from the methods of 
optimal control theory, in particular of Ref. 27, Chapter 3. 

Consider the case of equality constraints first. Clearly, the dimension 
of the control vector u must be greater or equal than 2 for the problem to be 
interesting, since for m = 1 the constraint (3.1) determines u(x,y) completely 
and there is no optimization problem left. 

We adjoin (3. 1) to the variational Hamiltonian with a set of Lagrange 
multipliers v(x, y) = ( v (x, y) , . . . , v (x, y)) as follows: 

T T 

H = L(z,u;x,y) + A (x,y)X(z,u;x,y) + (j, (x,y)Y(z,u;x,y) 

r- 1 ^ ^ 

+ v T (x,y)C(u;x,y) (3.4) 

The only change this brings about is in the optimality conditions: 
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3H 9L T T 9 ~ : 

— — = — — + X — + u, — + v 
9u 9u ~ 9u ~ 9u ~ 


9C 

9u 


(3.5) 


which, together with (3.1), represent p + r conditions to determine the p- 
component vector u(x,y) and the r-component vector v(x,y). 

r+J fv 

Now consider the inequality constraints (3.2). In the same way as above, 
we introduce the s-dimensional vector £(x,y) the components of which are the 
s Lagrange multipliers £ (x,y), . . . , £ (x,y), and define 

-L S 


T T T 

H = L + X X+ n Y + £ D 




(3.6) 


The necessary condition to be satisfied by this augmented Hamiltonian is: 
3H 9L T T .T 

9^r _ 9^r + ~ ^ + (3 - 7) 


which is the same as above with the additional requirement that: 


e^o, V o 

or: 

£ =0, D > 0 = 1, 2, . . . , s (3.8) 

Xj Xj 


Another approach is the following, again in the case of s inequality 
constraints: 


D (u;x,y )2 0 4 = l,2,...,s (3.2) 

Xj 


If we define H* as: 


T T 

H* = L + X X+ n, Y 
^ ,c ~ 


(which is nothing but the Hamiltonian defined for the case when there are no con- 
straints on the control variables), then, from (2. 6): 


6 J 


If IT 6 “ dxdy 

D ~ 


(3.9) 
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and this is, by definition. 


6J = ff 6 H*(z,u, A, ^ ;x,y)dxdy 


D 


where the A and the jj, satisfy the set of partial differential equations 


T T 

9A 9g 


9x 9y 


9H 

9z 


(2.7) 


with boundary conditions: 


,T T 9£ 

A P’ (or) a’(c r) = -~ 


(2. 7a) 


holding along the portions of 9D where z is not defined. 

For the control u(x,y) to be minimizing, we must have: 

<N/ 

6J s 0 

for all admissible 6u(x,y). This implies that we have: 

6H*sO 

for all admissible 6u(x,y) and this all over the domain D. 

Hence, at all points on D = 0 the optimal u has the property that: 


6 H* = - 6u s 0 

ou ~ 


(3.10) 


for all 6u such that: 

9D 

6 D = — — 6u 5 0, 
J} 3u ~ 


£ = 1, 2, . . . , s 


which yields (3.7). 

Another way of stating (3. 10) is to say that 5 H* must be non-improving 
over the set of all possible 6u. 
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Actually a much stronger statement, ”H* must be minimized over the 
set of all possible u " holds. This compact statement is the extension to two- 

r*-/ 

dimensional systems of the famous "Minimum (or Maximum, the distinction 

being due to a different convention in Russian literature and elsewhere to write 

28 

the Hamiltonian H*) Principle" . A rigorous proof of the above can be found 

29 30 

in a paper by Butkovskiy , in one by Egorov , and in the already referred 
25 

to paper by Lur' e , which presents an extension of the Weierstrass conditions 

of the classical calculus of variations to this type of problems and is nothing 

but an alternate way of stating the Minimum Principle. 

A remarkable approach using the function-space formulation is found in 
31 

a recent paper by Neustadt , who presents in the general non-linear case an 
abstract Minimum Principle covering almost all constrained minimization 
problems ever considered, including lumped and distributed variational problems 
under a variety of conditions. 

4 . Discontinuities in the Distributed Controls 

In the definition of the general problem given in section 1, we spoke of a 
"suitable class of functions". It is now appropriate to explain this notion, and to 
specify in detail the class of admissible controls and possible behavior of state 
variables. 

The distributed controls will be assumed to belong to a class of functions 

no less wide than the class of piecewise continuous functions in two independent 

variables. Possible discontinuities of distributed controls may occur along 

smooth, closed isolated curves A lying entirely inside the domain D. 

The state variables z. are assumed continuous across a curve such 

i 

that A ; the variables u^ are in general discontinuous , but their values on 
both sides of A are connected by the requirement that the tangential derivatives 
9z./ 8a along A be continuous, that is, if we denote by a subscript — the values 
of the variables attained on the curve A from one definite side of it, by + their 
values attained from the other side , 
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i 1,2, , 


(4. 1) 


[XT’ (tr) + Y.6 » (o-)]_ = [X.y' (cr) + Y.5 ’ (cr)] + , 

Recall that the curve A has the parametric representation: 
x = y(o) 

y = 6 (cr) 


n 


y and 6 being two piecewise continuous and differentiable functions of a. 

25 

It can also be shown that along such a line of discontinuity of distributed 
controls A: 


[\.S ' (cr) - y,.y' (cr)]_ = [A..6 ' (cr) - ^.y ' (cr)] 


i = 1,2, . . . ,n 


(4.2) 


[H 


^ 1 X_- tx T Y_] = [H-X T X_ - ^Yl 




/*w< 


(4.3) 


5. Remarks and Conclusion 

The necessary conditions have thus been established in a very general 
case of distributed-parameter systems. We are now confronted with a set of 
partial differential equations to be solved with adequate boundary conditions, 
which represent the first-order necessary conditions for an extremum. 

The original optimization problem is therefore reduced to that of the determi- 
nation of solutions for this resulting set of partial differential equations. This latter 
task unfortunately proves to be an extremely difficult one when we try to solve 
this system for some specific cases. The reason is clear: even for the 
simplest of all problems, the system formed by the necessary conditions will be 
greatly complicated. We start from one single partial differential equation as a 
constraint. This equation has to be broken, as shown in section 1, into a set of 
first-order equations, having now n variables z and p controls u.. After 
writing down the necessary conditions, we are left with a system of 3n+p 
partial differential equations in 3n+p unknowns. It is therefore understandable 
that only in very special cases can we find a. rigorous analytical solution to a 
given optimization problem. One such problem of structural interest is presented 


2k 



and solved in the next part of this study. The approach to the more complex, 
but at the same time, more realistic problems presented in parts C and D 
seems to have to be numerical at this stage of their solution. 
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PART B 


MINIMUM-MASS DESIGN OF A SIMPLY-SUPPORTED SHEAR PLATE 
FOR A FIXED NATURAL FREQUENCY 


A plate-like structure with special properties will be described and 
analysed. 

Given the shape of a homogeneous "shear plate" with constant thickness, 
we find its fundamental frequency of vibration. The problem considered is then 
to determine the thickness to be distributed within the same boundary so as to 
keep the above fundamental frequency constant while achieving the minimum 
possible mass. Only for simply supported edges will a complete solution be 
given. 

Stating the problem in mathematical form, we are confronted with a 
situation of the kind investigated in part A, for which an exact analytical solution 
will be found in the case of a rectangular or a circular boundary. 

1. The Shear Plate: Definition and Analysis 

We define a shear plate as a rather artificial plate-like structure such 
that the bending rigidity for normal loads is negligible: the running load is thus 
borne by shear only. 

We refer the plate to a set of rectangular Cartesian axes Ox, Oy in the 
middle plane of the plate; the plate occupies a simply connected domain D with 
a smooth boundary 9D. 

Now let us consider a rectangular element with sides parallel to the 
coordinate axes and with lengths dx and dy respectively, around a point 
P(x,y). Let h(x,y) be the thickness of the plate, p(x, y) its density, both at 
P, and let w(x,y;t) be the normal deflection of P at time t. We choose the 
axis Oz directly perpendicular to Ox and Oy. 

The forces acting on the above element (see figure 1) are all parallel to 
the z-axis. They are, in projection on Oz: 
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a) the shear forces: 


8(T xz h) 

-t hdy and + 1 r h + — dx J dy 

xz J \ xz dx ' 


on the two sides perpendicular to Ox. 


9(T vz h) 

-t hdx and + I t h + — t- 

yz v yz ay 


dy) dx 


on the two sides perpendicular to Oy. 
b) the d'Alembert force: 


,9 w , , 
-ph — — dxdy 

d€ 


From equilibrium, 


9(hr ) 8(hr ) 2 

xz yz , 9 w 

+ - p h~— =0 


dx 


9y 


8t 


But the classical stress-strain relations read 


T = Gy = G 
XZ XZ dx 


x = Gy = G 
yz yz 8y 


so that the differential equation of motion takes the form: 

2- (h ^ ) + _a. h ^ = o 

8x dx dy'dy G ^ 2 


( 1 . 1 ) 


This equation must be satisfied inside the two-dimensional domain D 
bounded by the frontier 9D. Let us now examine the boundary conditions. 

If the plate is supported at points (x^, y^) on segments of the boundary 
9D, we have the boundary condition: 
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w = 0 at (X^y^ 


(1.2a) 


If, on the other hand, the plate is free to deflect transversely at some 
other points (x , y ) on different segments of the boundary, there cannot be any 

Ci Ci 

shear component in the transverse direction. Consequently the boundary con- 
dition at such points is: 

G S = ° at 'v y 2> 

or 

fr = 0 at (x 2 ,y 2 ) (1.2b) 

The boundary curve 3D contains all points such as (x , y ) and (x , y ) , 

1 1 Ci Ci 

and only one of the boundary conditions (1.2a) or (1. 2b) must be satisfied at a 
given point on the boundary. 

We will confine ourselves to the case of a homogeneous plate (p = 
constant). Let us investigate the nature of the free vibrations of the shear plate 
with constant thickness. 

With h = constant, the equation of motion (1. 1) becomes: 

GV 2 w = p (1.3) 

ar 

which is in the same form that the equation of the free vibration of a membrane 

where G plays the role of the uniform tension T of the membrane. 

32 

Following the classical procedure , we write the displacement w as the 
product of a function W of the spatial coordinates only and a function f depend- 
ing solely on time: 

w(x,y;t) = W(x,y) f(t) 

(1.3) becomes: 

G \?W _ f 
p W f 


29 


where the superscript (') denotes a derivative with respect to time, and the 

2 

common value of these two ratios has to be a constant, -<o , negative for 
reasons of stability, f is then harmonic with frequency to , and W satisfies 
the partial differential equation: 

-GV 2 W = co 2 pW (1.4) 

inside the domain D, together with the boundary conditions: 


W = 0 

at 

(Xj.yj) 

(1. 5a) 

|W =0 

9n 

at 

( X 2’ y 2 ) 

(1. 5b) 


The eigenvalue problem represented by equation (1.4) together with the 
boundary conditions (1. 5) is of the classical type: 

L[W] = 1M[W] 

where the operators L and M are respectively: 

L = -GV 2 

M = p 

with boundary conditions independent of the eigenvalue X. The problem is easily 
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shown to be self-adjoint : if W and W are two eigenfunctions corresponding 

r s 

to two distinct eigenvalues X / X , then they satisfy the orthogonality relation: 

r s 

J'J’ pW^W^dD = 0 , r^s. 

D 

Also, it can be shown, as in the classical membrane case, that if the 
boundary condition (1. 5a) is satisfied over part of the curve 9D, the system is 
positive definite, and that if the boundary condition (1. 5b) is satisfied over the 
entire length of 9D, the system is only positive. 

To pursue the discussion further, we shall confine ourselves to an 
investigation of the properties of rectangular and circular plates only. 


30 




1. 1 Vibra tions of Homogeneous Rectangular Shear Plate with 
Constant Thickness 

The plate extends over a domain D defined by Osxsa, 0 s y s b. The 
boundary 9D is made of the straight lines x=0,a and y=0,b. Let: 



and write (1.4) as: 

V^w + a 2 w = 0 (1.6) 

where the expression for the Laplacian in rectangular coordinates is: 



(1. 6) is solved by separation of variables; to this end, we let the solution 
have the form: 

W(x, y) = X(x) Y(y) 

and, after substitution in (1.6) and rearrangement, obtain: 


_1_1M + _l_ d I M + 2 = 

X(X) dx 2 Y(y) dy 2 


which leads to the equations: 


^Vx<*> = 0 

dx 


d Y jp + y 2 Y(y) = o 
dy 

where 


2 2 2 
P +y =ol 
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so that, introducing four constants of integration A, A, , the solution W 

1 2 3 4 

has the form: 

W(x, y) = A sin |3x sin yy + A sin |3x cos yy + A cos px sin yy 

1 £j O 

+A 4 cos { 3 x cos yy (1. 7) 

where A^,A ,A ,A as well as p and y must be determined by means of the 

I Z O 4 

boundary conditions. 

a) Plate Simply-Supported Along the Edges 
The boundary conditions then read: 

W(0,y) = W(a,y) = W(x,0) = W(x,b) = 0 (1. 8) 

The first and third one of the above conditions lead to: 


A 2 = A 3 = A 4=° 


whereas the second and the fourth imply: 
sin pa = 0 
sin yb = 0 

and yield an infinite sequence of discrete roots: 


mu 

P = . 

m a 


m = 1,2, .. . 
n = 1,2,. . . 

so that we are led to the eigenvalues solution of the problem: 


rnr 

Y n = T’ 


a 


mn 




IT 




The corresponding modes are: 


m,n = 1,2,... 

(1.9) 
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( 1 - 10 ) 


. . mux rnry 

W = A sm sin 

nan nan a a 

The fundamental frequency of vibration is obtained for m = n = 1, and has 
the value: 


GJ — (V 

11 “ll 


Jf-M ¥9 


(1. H) 


b) Plate Simply-Supported Along Two Parallel Edges and Free Along the 
Other Two 

The boundary conditions now read, assuming the plate to be simply- 
supported along the parallel edges x = 0 and x = a: 


W(0,y) = W(a,y) = 0 


3W m 9W. .. n 

ay(x,o) = — (x, b >»o 


( 1 . 12 ) 


The integration constants in (1,7) are determined then from the first 
boundary condition: 

a 3 .a 4 .° 

From the third one: 


A 1 = 0 


and from the two remaining ones: 
sin pa = 0 
sin yb = 0 

P and y have to take the discrete values: 



mr 

b 


m = 1,2, 

n = 1,2, . . . 
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The frequencies of vibration are the same as in case a) above. The 


modes are now being given by: 


W = B 
mn mn 


. mirx niTy 

sin cos — — , 

a b 


m,n = 1, 2, . . . (1. 13) 


c) Free-Free Plate 

For this idealized case (realizable in practice by assuming the plate to 
be simply-supported at its center) , the boundary conditions read: 


8W. 3W _ . 
97 (x, 0 )=9^(x,b) 


0 


aw aw 

ir ( 0’y) = i7(a*y) = 0 


( 1 . 14 ) 


The natural frequencies of vibration are the same as in case a) or b). The 
natural modes are now given by: 


__ _ _ mirx mry 

W = C cos cos — — 

mn mn a b 


m,n = 1,2, . . 


(1. 15) 


1. 2 Vibrations of a Homogeneous Circular Shear Plate with Constant 
Thickness 

The circular plate extends over a domain D defined by 0 s r s a; the 
boundary of the domain is the circle BD described in polar coordinates by the 
equation r = a. 

With polar coordinates r and 0, the differential equation to be satisfied 
throughout D is: 

V 2 W(r, e) + o: 2 W(r, 0) = 0 (1. 16) 

where 

„ 2 
2 up 
ci = — 

“ G 

Assuming a solution of the form: 
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W(r, 9) = R(r) © (0) 


and recalling that in polar coordinates the Laplacian is given by: 


2 _ 8 2 ] 1 8 | 1 8 2 
8r 2 r 8r r 2 80 2 


(1. 16) reduces to: 
,2 


d“R 1 dR \ d 2 e ^ 2^ 

— + © +~— +aR© =0 

dr r dr / r d9 


which can be separated into two equations: 


d 2 0 2 

— — +m 0=0 


de 


d 2 R 

dr 2 


1 dR , / 2 m 2 \ 

V* ~ 2 ) 

r dr * r * 


R= 0 


(1.17) 


(1. 18) 


2 

The constant m is chosen to obtain a harmonic equation for © ; further- 
more, the solution must be periodic in e with the period 2ir; m has then to be 
an integer (m = 0, 1,2, . . .). 

(1. 17) yields as solution: 

© (9) = C„ sin m0 + C n cos me, m= 0,1,2,.. 

1 tyi Vm 


'lm 2m 

(1. 18) is a Bessel equation and its solution is: 
R (r) = C Q J (ar) + C, Y (ar), 

rn ^ ' 3m rr> ' ' 4m m ' ' 


m = 0,1,2,... 


where J (x) and Y (x) are the Bessel functions of order m of the first and 
m m' 

second kind, respectively. 

The general solution of (1. 16) can then be written: 

W (r, 0) = A, J (ar) sin me +A_, J (ar) cos m0 +A C Y (ar) sin m6 
m ' lm m v ' 2m m' ' 3m m 


+A, Y (ar) cos me , 
4m m 


m = 0, 1,2, . . . 


( 1 . 19 ) 
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a) Simply-Supported Plate 
The boundary condition reads: 

W(a, 0) = 0 

At every interior point of the plate the displacement must be finite: but 
Bessel functions of the second kind tend to infinity as the argument approaches 
zero. It follows that: 


A„ = A . =0 

3m 4m 

so that (1. 19) reduces to: 

W (r,0) = A„ J far) sin m.0 +A„ J (ar) cos m0, m = 0,1,2,... 
m lm m 2m m 

At r = a we must have: 


W (a,e) = 0 , 
m 


m = 0, 1,2, . . 


regardless of the values of 0, which implies: 

J faa) = 0 , m= 0,1,2,. 

m 


This equation represents an infinite set of characteristic or frequency 

equations, and for each m there is an infinite number of discrete solutions 

a . The fundamental mode of vibration is obtained for m = 0; as the first 
mn 

zero of Jq(x) occurs at x= 2.4048, the fundamental frequency is given by: 


u = 2. 4048 


G 


pa 


and the corresponding mode by: 


W 01 = A 01 J 0 < 2 - 4048 I> 


( 1 . 20 ) 


( 1 . 21 ) 


For the upper frequencies <o there are two corresponding modes to 
one particular frequency: the modes are then degenerate. 
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b) Plate Free Along the Edge 


The boundary condition now reads: 

W(r,0)| =0 

* r=a 

As above, the displacement at the center has to be finite, so that the 
solution is in the form: 

W (r, 9) = A, J (ar) sin m0 + A„ J (ar) cos m0, m= 0,1,2,... 
m lm m 2m m 

The boundary condition thus requires: 

J 

— J (err) I = 0 
dr m l r=a 

34 

But, by a property of Bessel functions : 


^ J 0 (ar) = - 0dJ 1 (o;r) 




m = 1,2, . . . 

the first formula being only a particular case of the second one if we note that: 


J (x) = (-l) m j (x), 
-m m 


Thus we must have: 


m = 0, 1, 2, . . 


J 

m- 


1 


(oa) - 


m 

aa 


J (aa) = 0 
m 


This equation represents an infinite set of characteristic equations, and 

for each m there is an infinite number of discrete solutions a corresponding 

mn 

to the zeros of this transcendental equation. The lowest frequency corresponds 
to m = 0 and is obtained from the first zero of J^(x), occuring at x = 3. 8317; 
thus, 


« ol = 3.8317 


G 


pa 


and the corresponding mode is: 


( 1 . 22 ) 
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(1.23) 


V’^^oiV' 8317 ^ 


Again, for each frequency w , m^O, there are two different modes, 


mn 


and it follows that for m ^ 0 the natural modes are degenerate. 


After this preliminary study of the vibrations of a shear plate, we now 
are ready to set up the optimization problem. We note that in the case of a 
circular boundary the optimal thickness distribution will evidently be rotation- 
invariant (symmetry about 0 and about any axis through 0), i.e. , independent of 
0, and the optimization problem will then reduce to the classical case of one 
independent variable, in this case r, encountered in optimal control theory. 


2. Optimization of a Simply-Supported Shear Plate 

Let us consider a simply-supported shear plate extending over a domain 
D bounded by the curve 3D. We assume that the material constituting the plate 
is homogeneous with density p and that the thickness is a constant h Q ; the mass 
of this reference plate is thus: 

M o = P h o Q 

where G is the area enclosed inside the boundary 3D. The fundamental 
frequency of the plate is as found above in the cases of a rectangular plate 
and of a circular one. 

We then want to find the optimal thickness distribution h(x,y) relative 
to an orthonormal set of coordinates such as to yield a minimum of the surface 
integral: 

M = p JJ" h(x , y)dxdy 
D 

or, equivalently, the surface integral: 
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( 2 . 1 ) 


J = 


If h(x, y)dxdy 

D 


the constraint the first frequency of vibration be 03 being expressed under the 
form of the partial differential equation: 


8 8w \ , 8 /U 9W , £L 2 U A 

ax v ax ay ay g f 

to be satisfied inside D, together with the boundary condition: 
w(x,y) = 0 on 3D 


( 2 . 2 ) 


(2.3) 


2. 1 Necessary Conditions for an Optimal 

Following the method derived in Part A , we break the partial differential 
equation (2.2) into a system of first order partial differential equations, with the 
help of auxiliary variables z. and u.. 

Let: 


Z l = 


Z 2 = 


Z 3 


w 


3w 

ax 

, dw 

h - — 

9y 


The system (A. 1. 1) then reads: 

a h_ 

3x h 
9z. 


9y 

^2 

9x 

Qy 


_3 

h 


= u. 


(continued) 
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(2.4) 


(2.4a) 


(2.5) 


( 2 . 6 ) 



(2.6a) 


(<r) ~ V-2 a ' ( CT ) = 0 
^ 3 P' (<*) - U^a' (o') = 0 

on 9D, uniquely determined by the parametric representation: 
x = a(o-) 

y = P(o') 

and they reduce to, considering (2.6): 

Z 1 = o 

*2=^ = 0 (2.6a) 

on 3D. 

Now the first three equations in the system (2.6) are very similar to the 
first two and the last of the system (2.4). We are thus led to the assumption: 




where a is any constant. 

Note that this assumption is compatible with the boundary conditions 



on 3D. 

The control equation reduces to: 
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1 - 


2 2 2 

Z _2_ _f3_ +£ 2fl _ 

2 2 G “f q 

an an 


or, expressed in terms of the old variable w: 


,9w. z .Sw, P 2 2 

+ V + G “f” 


( 2 . 8 ) 


which is a nonlinear first order partial differential equation to be satisfied by w, 
together with the boundary condition: 


w = 0 


on 9 D 


(2.8a) 


(2. 8) is a necessary condition for an optimum; we will show in Part C 
that it is also sufficient, as it is nothing but the expression in the particular case 
of a shear plate of an extremely general sufficient condition applicable to a broad 
class of structural optimization problems. Thus the whole optimization problem 
reduces to solving this partial differential equation: if a solution to (2. 8) can be found 
that satisfies the boundary condition (2. 8a), then the solution of the original 
problem can be immediately derived from the knowledge of w. 

Equation (2. 8) shows us + that a has to be positive: let 
2 

a = c 


Also, let: 


R 2 i 2 

gj = k 

G f 


and write (2. 8) in parametric form, introducing an unknown function 0 (x,y), as: 

9w r 2 " 2 2‘ 

— = ye + k w cos 0(x,y) 

9w ./"£ TFT' . A/ 

— - = Vc + k w sin 9(x,y) 

9y 

From the fact that the value of the second mixed derivative of w is 
independent of the order of the derivations: 

+ w attains the value 0 on the boundary, which is part of D. 

42 



— k 2 w ~ cos e - Vc 2 +k 2 w 2 ^ sin e 

“2~li - 2 A * 

c +k w 


2,2 2 
c +k w 


, 2 9 w . „ 2,2 2 ae 

k w — sin 0 + v c +k w — cos 0 
9x v dx 


90 . 90 

cos 0 — + sin 0 — = 0 
9x 9x 

which is a linear first-order partial differential equation that has to be satisfied 
by 0(x,y). 
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The characteristic lines are given by the system : 

dx dy _ dQ 

cos© sin0 0 

and are therefore defined by 


-x sin© + y cos0 = constant 


The characteristics are therefore straight lines, along which 0 keeps a 
constant value. 

On the boundary 9D, w has the constant value 0; therefore, 

9w , 9w , 
dw s — dx + - — dy = 0 
dx dy J 


cos0 dx + sin0 dy = 0 


which shows that 0(x,y) has a geometric significance: it is the angle that the 
normal to the contour 9D, oriented outwards, makes with a direction parallel to 
the x axis. The characteristics, which are, as we recall, straight lines along 
which 0 keeps a constant value, are therefore the normals to the contour. At any 


point P interior to the domain, the value of 0 is that of the angle the normal to 
the contour drawn from that point P makes with the direction Ox (Figure 2). 
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0 


x 


Figure 2. Solution of the partial differential equation (2. 8) for a 
general domain D: definition and interpretation of the 

different quantities introduced. 
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At P, 


dw = 


. / 2 ,2 2 ’ 

V c + k w 

. / 2 .2 2 ' 
Vc + k w 


(cos© dx + sin0 dy) 
dp 


if p is the distance from P to the boundary 3D, measured along the normal 
to the boundary pointing outwards. 

Thus: 


dp = 


dw 


./ 2.2 2 ' 
V c + k w 

or, integrating: 


p - p 0 = ^ log(kw + \/c 2 + k 2 w 2 ' ) , 


which is: 

w = ~~ sinh[k(p-p o )], e = ± 1 

For a point on the boundary (p = 0) , then w = 0: this shows that the con- 
stant Pq has the value zero, and the solution to (2. 8) is simply: 

w = ^ sinh(kp) (2. 9) 

where we recall that p is the distance from P(coordinates x and y) to the 
contour 3D. 

For a given contour, there are usually many normals that can be drawn 
from a point interior to the domain it encloses (figure 2). Thus the question 
arises to decide which normal to consider to define the distance p: it seems 
logical (and turns out to be the only solution for a rectangular plate) to take for 
p the shortest distance from P to the boundary 3D. 

Another approach to equation (2. 8) that reduces it to a form encountered 
in the theory of plastic torsion of a cylindrical beam is given in Appendix 2. 


**5 


We now apply the above result, valid for any domain D under the general 
assumptions outlined at the beginning, to a rectangular, and then to a circular 
shear plate. 

From a point inside a rectangular domain, we may draw four perpendiculars 
to the boundary. Therefore we have to consider 4 domains labelled (1), (2), (3), (4) 
as represented in Figure 3, delimited by the boundary, the bisectors of the 4 
right angles at the corners and the segment joining their points of intersection 
two by two (the latter reducing to a point in the case of a square plate). In each 
of those domains, the shortest distance of a point to the boundary is readily 
expressed, and this leads to an analytical expression for w, taking a different 
form in each domain*. 

Let the boundary be formed of the lines x = 0, x=a, y = 0, y = b. We 
may always assume, without any loss of generality, that a & b. Due to the 
symmetry of the problem, we need only consider one-fourth of the plate, for 
instance the portion delimited by the lines x = 0, x = a/2, y = 0, y = b/2. 

In region (1), w is a function of x only: 

w = ^-7 sinh(kx) 

and the optimal thickness distribution is found from (2. 2) which reduces to: 

0 

— | [h cosh(kx)] + kh sinh(kx) = 0 

The general solution to this partial differential equation of a special type 
is readily found to be: 

w (y) 

h= 

cosh (kx) 

*Note that it is easy to see why the shortest distance to the boundary has to be 
chosen in this case; another choice for p would not allow either the condition 
w = 0 to hold along the boundary 8D or the displacement w to be continuous 
inside the domain D. 
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Figure 3. The division of the rectangular plate into four regions, and 

the corresponding expressions for the optimal displacement w 



where «^(y) is an arbitrary function of y. 

Similarly, the optimal thickness distribution in (2) is given by: 

“ 9 <x) 

h= i 

cosh (ky) 

Now the shear force has to be continuous when we cross the boundary from 
region (1) to region (2); this means that the expression: 

, 9w . =* -> 

h — = h Vw • n 


has to be continuous across any line with normal n. In particular, between (1) 
and (2), we must have: 


, ,9w 9w. 

h( 9x " 9y * 


( 1 ) 


9w 9w. 

h( SF - i7> 


( 2 ) 


or: 


9w . 


9w | 

311 Id) 


9y | 

But: 




9w 


9w 


9x 

(1) 

3y 

1(2) 


(2) 


as the boundary line is defined by the equation y = x. 
Thus it is necessary that: 

h = 0 

on the segment y = x, 0 ^ x ^ b/2, which implies that 



and 

h s 0 
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The optimal solution is therefore the trivial one, zero. This result 

should not surprise us , as we already encountered it in one-dimensional 
7 13 14 

cases ’ ’ The trivial solution is a consequence of the constraint equation 

14 

being homogeneous in h. We also know from past experience that introduction 
of a minimum thickness constraint will not be of any help. 


2. 2 The Structural Mass Hypothesis 

The way to overcome this difficulty is to assume that the total mass of the 

plate is made up of two parts with drastically different structural properties: a 

constant fraction 6 is non-structural, whereas the remaining part of the mass, 
& 

labelled structural, is originally in the proportion 5 . 

The thickness is expressed as: 

b(x,y) = 5 1 h*(x,y) + 8 2 

where 

W l 

h* will often be referred to as the "virtual thickness", by opposition to 
the actual thickness h. 

Note that the non-structural mass is not at the control of the designer, 
and that we are trying to minimize the total mass by acting on the intermediate 
variable h*. 

The problem is then set up as follows: m inimi ze the total mass 


M 


= p fj h(x,y)dD = p6 x JJ h*(x,y)dD + p6 2 G 


D 


D 


or, equivalently, the surface integral: 


J* = ff h ’ lx , y)dxdy 
D 


( 2 . 10 ) 



with the constraint: 


— (h* — ) + — (h* — ) + £ oy (8 h* + 6jw = 0 
8x y 8x dy K 8y' G f v 1 2 1 


w = 0 on 3D (2. 11) 

The optimization process is similar to the one already performed; under 

exactly the same assumptions on A , A. , {j, , p, we find that the optimal displace- 

1 Z 1 o 

ment mode w has to satisfy the first-order partial differential equation: 


,3w 2 3w 2 2 p r 2 2 

'3x ; '8y ; Gif 

together with the boimdary condition: 

w = 0 on 3D 


( 2 . 12 ) 


which is exactly the same equation as (2. 8), with k now having the value 


P c 2 
Gif 


The solution is then, using the same notations as before: 


w = ~~ sinh(kp) 


(2.9) 


For the rectangular plate, the 4 regions previously defined (figure 3) 
have to be considered again; the optimal displacement w in each of them is again 
given by the same expressions, k having the new value 


co„ 



6 


1 


In region (1), from (2. 11), the optimal "thickness” h* has to satisfy 
the differential equation: 

3 6 2 

— [h* cosh(kx)] + k(h* + — ) sinh(kx) = 0 (2. 13) 


or: 
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3h* 

9x 


+ 2k tanh( kx)h* + k — tanh( kx) = 0 
6 i 


the general solution of which is: 

5 2 “l (y) 

h*(x,y) = + 

26 ^ cosh (kx) 

where S^y) is an arbitrary function of y. 

Similarly, in region (2), the optimal "thickness" distribution is given by: 

5 P w (x) 

h*(x,y) = - — + 

26 ^ cosh (ky) 


Now, as before, h* — has to be continuous across the line y = x, which 

an 

leads to the condition 


h* = 0 


on the segment y = x, 0 x ^ b/2. 

The expression for h* in region (3) will be obtained from that in region 
(1) by changing x into a-x; for region (4), from that in (2), y being replaced 
by b-y. The continuity of the shear from (2) to (4) when we cross the segment 

y = b/2, b/2 £ x a-b/2 

necessitates: 


h* 


3w 

8y 




(4) 


or, as: 


9y (2) 


9 w 

9y 


(4) 


/ 0 


on that line, 

h* = 0 on y = b/2, b/2sx^a-b/2. 


(2. 14) 
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Thus h* has to be zero on the lines serving to define the 4 regions 
which are not boundary lines (dotted lines in fig. 3). 

Applying this condition, we find the expressions for the undetermined 
functions w and to ; due to the condition (2. 14), region (2) has to be subdivided 

X m 

into (2) and (2 1 ) by the lines x = b/2, x = a-b/2 as shown in fig. 4. In the same 
fashion, those lines divide (4) into (4) and (4'). 

The optimal "thickness" distribution h* is: 


in(l) 


h*(x,y) = 


26 . 


cosh (ky) 
2 

cosh (kx) 


- 1 


in (2) 


h*(x, y) = 


26 


cosh (kx) 
2 

cosh (ky) 


- 1 


in (2*) 


h*(x, y) = 


26 . 


r cosh 2 (k|) 

2 

cosh (ky) 


- 1 


(2. 15) 


The expression for h* in (3) is obtained from that in (1) by changing x 
into a-x; in (4) and (4 1 ) respectively by that in (2) and (2* ) by changing y into 
b-y. We note that h* is nowhere negative, as it should. 

This corresponds to the actual thickness distribution: 


in (1) 


h(x,y) = — 


cosh (ky) 

x 2 
cosh (kx) 


in (2 r 


h(x,y) = 


1 + cosh (kx) 
2 

cosh (ky) 


(2. 16) 


in (2') 


h(x,y) = 


1 + 


,2 b 

cosh (Is.-) 

2 

cosh (ky) . 


the expressions in regions (3), (4), (4*) being derived as above. 

The contour lines (lines of constant thickness) corresponding to the case 
where 40% of the mass is initially structural and thus allowed to vary are 
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represented in figs. 5 and 6 respectively for a square plate (a/b = 1) and for a 
rectangular plate for which the ratio of the larger side to the smaller one is 
equal to 1. 5 (a/b = 3/2). We note that the optimal thickness distribution attains 
its min imum of 0. 6 corresponding to h* = 0 along the lines already described 
in fig. 3, and that a thickness of more than 1. 7 times the original one is attained 
at the four points which are the middle of the sides in the case of the square plate. 

The optimal mass of the plate is given by the double integral: 

M = p //«* , y)dxdy . 

D 

We compute it for one-fourth of the plate: due to the symmetry, the total 
mass will be four times that result. 

We obtain, after performing the double integration, the exact expression: 

, , c r ab 2 . , 2 „ b , a-b . , _ , . , 

M = P 6 2 V~7T + ~2 sinh < k 2 ) + 2k" smh<kb)] 
k 

where: 


k = < °£\/f^ = ' r ^ l (a2+b2) 

This leads to an optimal mass ratio equal to: 

2ab 


to = < - + 


2 12 2. 2 2 
v v 6 ^(a +b ) 


a-b 




2ir \J 6 (a 2 +b 2 ) 


sinh 2 ~ \/s 1 (a 2 +b 2 )' 

i \ 1 (a 2 +b 2 )^j 


(2. 17) 


As could have been expected, this quantity only depends, regarding the 
geometry, on the ratio a/b = X of the lengths of the two sides of the rectangular 
plate, and fo is rewritten as: 



0 


a 


Figure 5 


a_ 

2 


Optimal configuration of a square simply- supported shear 
plate (the thickness of the reference plate is taken equal to unity). 
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8 , = 0.4 
S 2 = 0.6 



x 


Figure 6. Optimal configuration of a rectangular simply- supported 
shear plate (a/b =1. 5). 



iTl = 6 


=■ v/mAi! 


2 / V - 81,11,2 s V6 i <x - +1) 

6 1 (X + 1) 


0 + ~ 

— — ^ sinh t \Z& „ <A 2 + 1) * 

/ 2 ^ A 1 J 

2ir x/e^A +1 > 


(2. 17) 


Let us now examine the extreme cases: when 6 -• 0, i. e. , when all 

the mass is non-structural and thus impossible to redesign, the optimal mass 
ratio tends to the value 1, as it should, showing that the result of the optimization 
process is the original plate itself, as can be seen by noting that in the vicinity 
of 6^0, 

2 „ ,.2 ... 

— ; V 77 6 a +1) 

<* +D — 2 

4A 2 

sinh ^ \/5 1 (>- 2 + 1) ~ ^ \/6 1 (X 2 + 1)' 

On the other hand, for 6=1 and 6=0 (case examined first when all 
the mass is structural and allowed to vary), this ratio is equal to zero, due to the 
fact that the optimal thickness distribution is then identical to zero throughout the 
domain D. 


sinh 


— \/6 
2A v 1 


For a square plate, a = b thus A = 1, and the optimal mass ratio is then 
equal to: 



(2. 17a) 


For the shapes pictured in figs. 5 and 6, the mass savings are found equal 
to 14. 3% and 15. 9% respectively, which is an encouraging result if we consider 
that 60% of the initial mass is not at the control of the designer. 

When A goes to infinity, i. e. , when one dimension of the plate becomes 
extremely large compared to the other one, then the optimal mass ratio tends 
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to the limit value: 



which is identical with the optimal mass ratio found for the minimum-mass design 

14 

of the one-dimensional cantilever wing for fixed torsional frequency . The 
analogy between that problem and our shear plate one is by the way worth noting. 

Also, for a fixed value of 6 , the mass ratio first decreases when the 
ratio a/b increases, reaches a minimum and increases again to the value to 

CO 

when 6^ goes to infinity. This variation of In with X for a 6^ of 0. 3 is plotted 
in fig. 7. The asymptote is the straight line to= 0.901663, and the maximum 
mass saving is obtained in that case for a/b = 3.42. The value of the mass ratio 
is then 0. 88621, corresponding to a saving of 11.4% as compared to a saving less 
than 9% for a square plate made of the same material. 

The variation of the mass ratio to with the proportion 6 ^ of structural 
mass is represented in fig. 8 for different values of the ratio a/b. 

2.3 Application to the Circular Plate 

We now turn our attention to the optimization of a circular simply-supported 
shear plate. The optimal thickness distribution being assumed to be rotation- 
invariant, i. e. , its expression in polar coordinates to be independent of 0 , the 
problem might also, as pointed out in the introduction, be viewed as a one- 
dimensional one, and solved using the classical methods of optimal control theory, 
which will provide two ways of finding the solution. 

*f* 

In polar coordinates (r,9), an axisymmetric optimal solution h* has to 
satisfy (2. 2) rewritten as: 


where the superscript (*) has the same significance than before, the actual thick- 
ness being h = 6 h* + 6 , with 6 +6=1. 

1 ^ A. Ct 
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OPTIMAL MASS F 




OPTIMAL MASS RATIO 



8 , 


Figure 8. Variation of the optimal mass ratio with 6^ for different 
values of the ratio a/b. 


60 


(2.18) 


1 2 2 
(h*w' )' + — h*w' + k (h* + — — ) w = 0 

where (') denotes the derivative taken with respect to r. 

With the value of the optimal displacement given by (2. 9) where p is 
taken equal to a-r. 


w = — sinh[k(a-r)] (2. 19) 

the optimal virtual thickness h* has to satisfy the ordinary differential equation: 

h*' + (— - 2k tanh[k(a-r)]) h* - k tanh[k(a-r)] = 0 (2. 20) 


For the same reasons than for the rectangular plate (continuity of the 
shear), theoptimal virtual thickness h* has to vanish at the center of the plate. 
The solution is thus found to be, after some manipulation: 


2 sinh(2ka) - sinh[2k(a-r)1 - 2kr cosh[2k(a-r)] 

5 1 8kr cosh 2 [k(a-r)] 


( 2 . 21 ) 


and the optimal thickness distribution is given by: 


h= M r 


(5 


4 cosh [k(a-r)J 


sinh(2ka) - sinh[2k(a-r)] 
2 

8kr cosh [k(a-r)] 


) 


( 2 . 22 ) 


We recall that: 

k -“ f yis 


2.4048 



The optimal mass ratio is then equal to: 


a 

f 2Trhrdr 



Tra 
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and, after some easy integrations, found to be: 


in = 



sinh (ka) 

, 2 2 
k a 


) 


(2.23) 


As expected, the ratio is equal to zero when all the mass is structural 
(6 =1), and tends to the value 1 when all the mass is non- structural (6^ = 0), 

and therefore not at the control of the designer. 

A diametral cross-section of the optimal plate corresponding to the case 
where 60% of the mass is non-structural is represented in fig. 9. The mass 
saving then obtained is equal to 8. 6%. The variation of the optimal mass ratio 
with the proportion 6 of structural mass in the reference structure is plotted 
in fig. 10: note that the savings obtained for a circular plate are slightly less 
important than those obtained in the case of a rectangular plate, as can be seen 
by comparison of the curves represented in fig. 7 to the one of this figure. 

To solve the same problem as a direct application of the methods of 
optimal control theory, we have to minimize the mass: 


M 




hrdr 


or, which is equivalent, the definite integral: 


a 


J = J" hrdr 
0 


(2. 24) 


The constraint is given by equation (2. 11) rewritten for the case of axi- 
symmetric vibrations, using polar coordinates , as: 


-| n 

(h*w' )' + — h*w' + *7w.(5h*+5)w=0 
r G f 1 2 


(2.25) 


(' ) denoting as before the derivative with respect to r, together with the boundary 
condition: 
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I 


Figure 9. Diametral cross-section of an optimal circular simply- 
supported shear plate for which 40% of the mass is 
structural. 
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w(a) = 0 


We introduce the auxiliary variable: 
s = h*r 

and rewrite (2. 25) as: 

2 6 2 

(sw 1 )'+ k (s + ~ — r) w = 0 


w(a) = 0 


where: 

.2 o 2. .2. 4048, 2 

k = 6 , = ( 6, 

G f 1 a 1 


Introducing another auxiliary variable t by the means of: 
t = sw' , 

we rewrite the optimization problem as: 

Minimize: 


/ 


J = J sdr 
0 

subject to the differential constraints: 
w’=^ 


5 

t’ = -k 2 (s + zr— r)w 
1 


and the boundary conditions 


(2.26) 


(2. 24) 
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w(a) = 0 


t(0) = 0 + 

The Hamiltonian is formed: 

5 

H = s + X — - k^X (s + ~ r)w 
w s t o 


and the necessary conditions for a minimum are given by: 


X' 

w 


X' 

t 


2 2 
k t (S + 6 r) 


X 

w 

s 


X w l 2 

! ” ~~2 ~ k \ W = 0 ’ 

s 


the last equation being the control equation. 

The natural boundary conditions are: 


X 

w 


( 0 )= 0 


X t (a) = 0 


A solution to the problem appears to be such that: 


w a 


X 

t 


w 

J 

a 


( 2 . 27 ) 


a. being a proportionality constant, assumption compatible with the original 


^The second of the conditions results from the fact that s(0) = 0 imposes the 
vanishing of t for r = 0. 
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boundary conditions. The control equation is then rewritten as: 


t , 2 2 n 
a - -£■ + k w =o 

s 

or: 

2 2 2 
w' = a + k w 

which shows that a has to be positive: let 


ol = c 

(2.28) is written as: 
dw 


dr = e 


J 2 ,2 2 

Vc +k w 


e = ± 1 


or, integrating: 


(2.28) 


e T „ , / 2 . 2 2\ 

r-r =— Log(kw + \Jc + k w ) 
ok v 

which is also: 

w = ^ sinh[k(r-r Q )] (2. 29) 

and, due to the condition w(a) = 0, the optimal displacement mode is finally 
given by: 

w = -^ sinh[k(a-r)] (2. 30) 

which is exactly the expression of the solution found to the two-dimensional 
problem, as rewritten for the case of axisymmetrical vibrations in polar 
coordinates (equation (2. 19)). The problem of finding the optimal thickness 
distribution from now on is the same as treated above. 

The solutions to the optimal problem obtained from two different methods 
— satisfying the necessary conditions for optimality in the classical one- 
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dimensional case and in the two-dimensional one, as derived in part A — are 
therefore the same when the two-dimensional domain can be in fact reduced 
to a one-dimensional one in an image-space as it is the case for the axisymmetric 
vibrations of a circular plate, therefore proving the validity of the optimality 
process in two-dimensions for this particular case. 


3. Simply-Supported Rectangular Shear Plate with a Minimum-Thickness 
Constraint 

As we have previously seen, the optimal thickness distribution in the case 
of a rectangular plate is found such that the virtual thickness h* vanishes along 
the segments constituted by portions of the bisectors of the four corner right 
angles and of the line joiningtheir points of intersection, as represented in fig. 6. 
This has as a consequence the instability of the structure, and a plate built 
following this theoretical distribution would immediately collapse if simply - 
supported. The same situation arises for a square plate, which is a particular 
case of the above for which h* is zero along the two diagonals. 

It is therefore very desirable to impose a constraint on the virtual thick- 
ness h*; imposing h* to be at least equal to h* is of course equivalent to 

impose a minimum actual thickness h = 6 h* + 6 . 

o 1 o 2 

The optimization problem is the same as before: 

Minimize the surface integral: 


J* = J'J' h*(x,y)dxdy 


(3.1) 


D 


subject to the constraint: 

-^-(h* "Ip-) + -~(h* — ) + -^w^(6 h* + 6 ) w = 0 
9x v dx dy dy G f' 1 2 

w = 0 on 8D 

with the additional inequality constraint on the control variable: 


(3.2) 
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(3.3) 


h* - h* SO 
o 

(3. 2) is broken down into a system of the form (A. 1.1). Following the 
method derived in (A. 3), we form the augmented Hamiltonian: 

z z 

H = h * + \il + Vl + X 3 U 3 + ^li^ + ^2 u 2 


-^ [u l + G“f 2(5 l h *- f6 2 )2 l 1+e<h *' h ? 


(3.4) 


where: 


£<x,y) < 0 
£(x,y) = 0 


if 

if 


h* = h* 
o 

h* s h* 
o 


(3.5) 


The necessary conditions for an extremal are expressed by the system of 
partial differential equations: 

8x h* 

9y h* 

dZ 9 

sr-"i 


dz r 


U, 


ay 2 

9x 3 


— — = - u u> 2 ( 6 h* + 6 )z 
8y 1 G r 1 2' 1 


(continued) 
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(3.6) 


(3. 6a) 


(continued) 
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a being a proportionality constant. The control equation reduces to, going back 
to the original notations: 


,3w.^ ,9w.^ o 2 C 2 t 

<S> + v = “ + c“f 5 l w +5<X ’ y) 


(3.7) 


At every point of the domain D where h* is greater than h*, 


£(x,y) = 0 

and (3.7) reduces to: 


(—) 
v 9x ’ 


+ (— > 2 

3y 


0 ! + 


Si 2 c 

GO 0 

G f 



Let 


2 

a: = c . 

The solution to this equation was found to be: 
w = sinh(k p) 


(3.8) 


where: 

k= “f JlFi 

and the quantity p having been previously defined. 

We now focus our attention onto the case of a rectangular plate with sides 
of lengths a and b (asb), For obvious reasons of symmetry, it is sufficient 
to consider the portion of the plate which is the left half of the region previously 
called (2). 

Inside this portion, p is equal to y, and: 

w = “■ sinh(ky) (3. 8a) 

K 
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The optimal virtual thickness is, as before: 


h*(x,y) 


6 

2 to(x) 

2 

26 ^ cosh (ky) 


(3.9) 


where To is an arbitrary function to be determined. 

Now in the regions where the minimum value permissible for h* is 
attained, that is: 


h* = h* 


then equation (3. 2) reduces to: 


2 

V w + 


R 2 
G f 



w = 0 . 


(3.10) 


This is an equation of the form (1. 6), with: 



(3.11) 


The general solution is, introduce the four constants of integration 



w(x,y) = A sin px sin yy + A sin px cos yy + A cos px sin yy 

-L Z (J 

+ A^ cos px cos yy (3. 12) 

where: 

2 2 2 

p +y =a . 

The problem is now to find the shape of the curve A dividing (2) into two 
regions with different properties: in the first one, which should be attached to the 
bisector of the lower-left corner and to the line defined by 

y = b/2 , b/2sxi a/2 , 

h* keeps the minimum permissible value h* ; in the second one, closer to the 
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side y = 0, the distribution of h* is as given by (3.9). 

The shape of this dividing curve A can be imagined quite easily, with the 

help of the shape of the lines of constant thickness already derived for the case 

where there is no constraint on the thickness, and from previous experience in 

14 

one-dimensional cases : A should consist of a curve starting off somewhere on 
the lower side of (2) between the points x = 0 and x = a/2, situated on the same 
side of the line y = x, and prolongated by a nearly straight line parallel to the 
x-axis after the point with abscissa b/2. In any case, in the corridor enclosing 
the bisector y = x, the displacement w has to be symmetrical about the same 
line, and this can happen if and only if: 



To determine the exact shape of A , we note that w and its first 
derivatives 9w/9x and 9w/9y have to be continuous across it. 

A point (| , q ) of A has to be such that: 



sinh(kq) = a sin p| sin pq + a (sin p| cos pq + cos p| sin pq) 
-L ^ 

+ a^ cos p| cos Pq 

0 = a. cos p| sin p-r) +a (cos p| cos p-q - sin p| sin Pq) 
1 ^ 

- a 4 sin p| cos Pq 


- cosh(kq) = a 1 


- a . 


sin p| cos Pq - a (sin p| sin pq - cos p| cos pq ) 
cos pi sin Pq 


(3.13) 


(the a , a , a replace the previous constants A , A , A as the displacement 

i. Z 4 JL Ct . t 

mode w is defined up to a constant multiplicative factor). 
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The point 0 belongs to the domain in which w is given by (3 . 12) . The 
condition w = 0 at the corner 0 leads to: 


= 0 . 


By subtraction of the last two equations in (3. 13), we obtain: 

k cosh(kr| ) 
a l (3 sin |3(|-r|) 

From the second equation, 

kcosh(k-n)cos p£ sin p-q 
a 2 p sin P(£ -t|)cos P(£ + q) 

Replacing a and a by their values in the first equation (3. 13), we 
find the relation between £ and r) : 


tanh(kri) = 


ksin(Pq ) 
Psinp(£ -r| ) 


[sin p£ - cos p£ tan P(£ + q)] . 


After some trigonometric manipulation, this equation turns out to be 
solvable in £ , yielding: 


1 . -1 F . , 2k sin (P"n ) I 

2p L 13 tanh < kr > ) J 

where sin ^ is defined by: 

I Arcsinu + 2hv 
S * n U j tt- Arcsinu + 2hiT , 


(3-14) 


h= 0, ± 1,. . . 


Arcsinu representing the principal determination of arcsinu, and where the proper 
choice of the integer h has been made. 

Recall that k and p are defined as: 


VFs - Sb v/ 6 i< a2+ b2 >' 


(continued) 
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(3 = o> 


/Z773Y . l A(T3) 

V 2G V 1 h*/ ab v 2 \ 1 h* ^ 


(a 2 + b 2 ) 


The curve A is thus defined by the equation: 

x = £ (y) • 

Although we cannot do the inversion in an analytical form, its Cartesian 
equation is, formally: 


y = n (x) = £ _1 (x) . 

We are now able to write down the expression for the thickness distribution 
in region (2): it has to be continuous across A; thus, at every point (£ ,p) of 
A: 


_ il. + = h * 

26 cosh 2 (kri) J 

where the unknown function w is defined as: 


O) (x) = 



cosh [kp(x)] 


On the right of and below A, the optimal "thickness" distribution is 
given by: 


h*(x,y) 




2 

cosh [kr| (x) ] 
2 

cosh (ky) 


(3. 15) 


(in region (2), the extension of this formula to the other regions being immediate), 
and this corresponds to a true thickness distribution of: 


h(x,y) 


6 


2 


2 


+ (6 x h* + 



2 

cosh [kp (x)] 
2 

cosh (ky) 


(3.16) 
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Of course, in the remaining part of domain (2), h keeps the constant 
minimum value: 

h(x,y) = 6^* + 6 ^ • 

The distribution in the other 3 regions follows immediately from 
symmetry. 

We cannot here give an explicit expression for the optimal thickness 
distribution, as the function ti(x) cannot be put under any analytical form and 
is only known implicitly. 

Of course a computational approach is always possible: contours of 
constant thickness have been represented in figures 11 and 12 respectively for the 
case of a square plate and that of a rectangular plate for which the ratio of the 
lengths of the sides is equal to 3/2, as in the previous unconstrained problem. 
Again, the initial proportion of structural mass in the reference structure 5 
was taken equal to 0.4, which means that 60% of the mass is non-structural. 

The constraint h* was taken equal to 0. 5, representing an actual constraint of 
0. 8 on the true thickness h. 

In the case of the square plate, 84.4% of the plate surface is at the 
minimum allowed thickness of 0. 8. The mass saving then obtained is of 14. 1% , 
slightly inferior to the one obtained in the case where there was no minimum 
constraint on the thickness. 

For the rectangular plate (a/b = 3/2), the minimum allowed thickness is 
reached on 77. 6% of the total surface. The total mass is found to be of 0. 845 
that of the original plate, corresponding to a saving of 15. 5% that compares well 
with the 15. 9% obtained when there was no constraint on the thickness. 

When we compare the corresponding optimal shapes for the cases where 
there is a minimal constraint on the thickness and the unconstrained ones, the 
main feature that we notice is the augmentation of the area of the domain formed 
by the reunion of the points where the thickness was inferior or equal to the 
given constraint, together with the drastic diminution of the maximal thickness, 
that was obtained at the four points middle of the sides in the case of the square 
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rolo- 



Figure 12. Optimal thickness distribution for a simply-supported 
rectangular shear plate with a minimum -thi ckn es s 
constraint (a/b = 3/2). 
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plate, and at the two points middle of the larger sides in the case of the rectan- 
gular shape. Physically, this is logical, as in the constrained case, the parts of 
the plate where the thickness is above the imposed minimum have less to con- 
tribute to the rigidity of the structure than in the unconstrained one, where they 
have to compete with the effect of antagonistic hinges. This can also be viewed 
in terms of the total potential energy for the first deflection mode of the plate, 
which has to remain the same for a fixed natural mode, whether or not there is 
a constraint on the thickness: the contribution of a large thickness at some points 
of the plate has to balance that of hinges in the unconstrained case, whereas 
everything gets smoother when a constraint on the thickness is imposed. 

This feature (diminution of the maximum thickness attained and augmenta- 
tion of the area where the thickness is inferior or equal to the given constraint 
as soon as it is applied) is the exact extension to two -dimens ions of properties 

already recognized for the one-dimensional problems where a minimum con- 

14 

straint was imposed on the thickness 

For the practical design, these results are of extreme importance. The 
shapes obtained when the thickness is subject to a minimal constraint are much 
easier to realize in practice: the maximum thickness attained in the domain 
covered by the plate is much smaller than in the unconstrained case, most of 
the surface of the plate is at this minimum imposed thickness, and there are 
no hinges as before. Moreover, together with all these advantages, the striking 
fact is that the mass saving is only slightly inferior than in the unconstrained 
case. 

14 

Again, these properties were noticed also for one-dimensional cases 
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PAST C 


THE MINIMUM-MASS DESIGN OF A SANDWICH PLATE FOR A GIVEN 
FUNDAMENTAL FREQUENCY OF VIBRATION 


A further step in the domain of applications to more realistic structures 
is that of the optimization of plate-like structures encountered in practical design 
for a given fundamental frequency of vibration. For a sandwich plate, the 
necessary conditions lead, when the displacement along the edge is assumed to 
vanish, to an optimality condition expressing the uniformity of the optimal 
energy distribution in the form of a second order nonlinear partial differential 
equation for the optimal displacement mode w that does not involve any design 
parameters. The optimal thickness distribution is the solution of a second 
order linear partial differential equation. 

A numerical solution of these equations is presented for a square aluminum 
alloy-aluminum honeycomb panel simply-supported along the edges. 

1. Statement of the Optimization Problem 

A sandwich plate is a composite plate consisting of a core layer of 
thickness h Q and of two face layers of thickness t. It is assumed that t is 
small compared to h Q and that the core material is much more flexible than the 
face material. Under these assumptions the transverse shears are predominantly 
taken by the core plate while the bending stresses are primarily taken by the face 
plate (Fig. 13). We will take the elastic modulus in transverse direction of the 
core-layer material, G , to be infinite. 

The core material is generally extremely light, and despite the small 
thickness of the face material, most of the weight is concentrated in the faces. 

We will therefore act on the face layers with thickness t, and our goal will be 
to find the optimal t distribution of the rectangular plate of minimum weight 
having the same fundamental frequency of vibration as a uniform reference one. 
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Under the small-deflection assumption, the equations of equilibrium for 

36 37 

a sandwich plate have been derived by Reissner ’ who followed the classical 

approach of combining the equilibrium equations and the stress-strain relations, 
38 

and by Hoff who used a variational approach. Assuming that the only external 
action applied to the plate is the normal force of intensity q, the equilibrium 
equation reads: 

a 2 M 9 2 m 9 2 m 

x „ xy y 

— — - 2 *- + — - -q (1. 1) 

9x dxdy 9y 


where M^, M are the bending moments per unit length of 
perpendicular to the x and y axes respectively, M the 
per unit length of sections perpendicular to the x axis , the 


according to the Timoshenko convention 


39 


From the stress-strain relations 


36 


sections of the plate 
twisting moment 
signs being taken 


M = - D 
x 



+ 


v 


9 2 w^ 

9y 2 ' 


M 

xy 


D(l-v) 


2 

9 w 
9x9y 


= - M 


yx 


M = - D 

y 




) 


where the bending stiffness factor D is defined as: 


( 1 . 2 ) 


t(h + t) 2 E 

D " ~ - 2 -~ 

2(l-v ) 

being the elastic modulus of isotropic face-layer material, 
ratio. 


(1.3) 

v the Poisson 


Substituting these expressions in the equation of equilibrium (1. 1), we 

obtain: 
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(1.4) 


V 2 (DV 2 w)+(1 


\ rTXrl' 


a 2 w 8 2 d ~ 2 


dxdy dxdy 9x 8y dy dx 


8 2 w 8 2 D 8_w\ 

a,, 2 "*, 2 a 2 / " q 


Under our hypothesis t « h^, a very good approximation for D is: 


D = 


2 

h E 
0 s 

2(1- v 2 ) 


(1.5) 


and (1.4) reduces to: 


V 2 (tV 2 w)+(l-v) 


\ Ffxi 


2 2 2 2 2 2 
' 8 w 8 t 8 w 8 


dxdy dxdy dx 2 dy 2 8y“ dx' 


t 8 w \ 

2 ^2 J 

Av / 


2d-_v 2 ) 

h 2 E 
0 s 


( 1 . 6 ) 


For a sandwich-plate of uniform face thickness t and uniform core 

thickness h , the differential equation for free vibration is (1. 6) where q is 

U 2 2 
replaced by the d' Alembert force -(h^p + t Q P s^ ) and t is a constant 

'0= 


4 >; = 2(l-v 2 ) h O P c +t O p s 8 2 w 

2 


V w = 


h 2 E 
0 s 


0 


£ 

8t 2 


(1.7) 


is the density of the core-layer material, p^ that of the face-layer 

material. 

Denote by n and s the coordinates in the directions normal and 
tangential to the boundary 8D. For a clamped edge, we have the geometric 
boundary conditions: 


3w 

w = 0 and — = 0 along 8D . 
8n 


(1.7a) 


At a simply-supported edge, the boundary conditions are: 


w = 0 and M = 0 
n 


along 8D 


(1.7b) 


where is the bending moment per unit length associated with the cross- 
section whose normal is n. 
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In the ease of a free edge, the boundary conditions are: 


3M 

ns 

M = 0 and V =Q - — =0 along 3D (1.7c) 

n n n 3s 


where V denotes the vertical force, Q_ is the shearing force and M is 
n n ns 

the twisting moment about the direction n. All three quantities V , Q and 

n n 

M are for one unit length of boundary and are associated with the cross-section 
ns 

whose normal is n. 

The relation between the moments and shearing forces and deformations , 

32 

in terms of normal and tangential coordinates , are : 


M 

n 


= -DV w + (l-v)D 


( 


3_ 3w 
R 9n 


+ 


3 2 w \ 
3s 2 ' 


M 

ns 


(l-v)D 


/ 

I 3 w 
\ 3n3s 


1_ 3w \ 
R 3s / 


Q 


n 



2 

w 


where the Laplacian has the form: 


V 


2 


3 2 + 1 _ 3 _ 

„ 2 R 3n 
3n 


3 


2 



and R denotes the radius of the boundary curve. 

To formulate the eigenvalue problem, we let in the classical fashion the 
displacement w be given by: 

w(x,y;t) = W(x,y)f(t) 

where W depends on the spatial coordinates only and f is a time-dependent 
harmonic function of frequency w. (1. 7) reduces to: 


h o E s 


Vc +t 0 p f 


2 

cj W 


(1.8) 
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_gndjhe^ boundary conditions are, depending on the situation at the edge, given by 

one of the conditions (1. 7a), (1. 7b) or (1. 7c) where w has been replaced by W. 

32 33 

This eigenvalue problem is a classical one ’ ; it is shown to be self- 

adjoint, and yields an infinite sequence of eigenvalues, the smaller one being the 
fundamental frequency of vibration of the plate under the imposed boundary 
conditions . 

The variable face-layer thickness t of a. rectangular plate having for one 
of its natural frequencies of vibration — that we can reasonably assume to 
be also its fundamental frequency if the plate does not differ too drastically from 
the uniform one — together with the corresponding mode w, have to satisfy 
the partial differential equation: 


V 2 (tV 2 w)+(l-v) 




rv2 J2 r\2 J2 r.2 

9 w 9 t 9 w 9 t 9 w 

3x3y 9x3y 3x 2 3y 2 3y 2 3x 2 


) 


11 tZ ) 

h o E 

0 s 


( Vc +t Pf )w f W 


= 0 . 


(1.9) 


If we introduce the non-dimensional face-layer thickness r = (t/t^), the 
above equation reduces to: 


V 2 (r V 2 w)+(1 


-v) ^ 


3 2 r 3 2 w 


3 2 r 3 2 w 3 2 t 8 2 w 


2 2 

3x3y 3x3y 3x 3y 


2 2 
3y 3x 


)- 


p (t + k) w = o 

( 1 . 10 ) 


where 'K and p are constants depending only on the material properties of the 
reference structure, defined as: 


k = 


*0 P[ 


( 1 . 11 ) 


2 2(1 - V) Pf 2 

P = X CO .. 

p 2 f 

h n E 
0 s 


( 1 . 12 ) 
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We are now able to state the optimization problem: mininize the total 


mass: 


M = ff (P c h Q + P f t)dD 


D 


or, equivalently, the functional: 


If 


J = II T dD 
D 


(1.13) 


taken over the domain D covered by the plate, subject to the partial differential 
equation constraint: 


0 


„2 _2 ( „ a 2 r 3 2 w a 2 r a 2 w ~ 2 ~ 2 

V (rV W)+(1-V) ( 2 — — 

9x9y 9x9y 9x 9y 9y“ 9x“ 


9 2 t 9 2 w > 

~ a 2 a / " 


p (t+k)w = o 


( 1 . 10 ) 

together with boundary conditions (1. 7a), (1. 7b) or (1. 7c) along the portions of the 
boundary 8D where the edge is respectively clamped, simply- supported or free. 

For a rectangular plate with sides of length a and b, a ^ b, the natural 
frequencies of vibration are found to be: 

1 


CO 


mn 


2 fm 2 n 5 

' (jr* + 


0 


2(l-v 2 ) Vc + V>f 


(1. 14) 


m,n = 1,2,3, .. . . 

The fundamental frequency is: 


CO = CO 

f 11 


=* 2 (V-) 

V 2 2/ 



2(1- v 2 ) Vc +t 0 Pf 


(1.15) 


and the value of (3 in equation (1. 10) is, in that particular case: 
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(1.16) 



The following derivation of the necessary conditions will be valid for the 
most general shape of the boundary 3D. For the applications, we will assume 
the plate to be rectangular and will take the coordinate axes to coincide with two 
of the sides, the origin being one of the corners of the plate. Only clamped and 
simply-supported edges will be considered. The boundary conditions will be 
picked among one of the following: 

i. edges parallel to the x axis. 

CLAMPED: w = 0, 3w/ 3y = 0 along y = 0, b 

2 

SIMPLY-SUPPORTED: w = 0, r = 0 along y = 0, b 

ay 2 

ii. edges parallel to the y axis: 

CLAMPED: w = 0, 9w/ 9x = 0 along x = 0, a 

2 

SIMPLY-SUPPORTED: w = 0, r = 0 along x = 0, a (1. 10a) 

3x 


2 • T he Necessary Conditions 

The constraint (1. 10) is transformed into a set of first order partial 

differential equations in the form (A. 1. 1) following step by step the general 

2 . 

method outlined in Part A. 9 r/9x9y is computed from (1. 10), leading to the 
system: 

z^ = w 


9z 


1 


9 x 


(=— ) 
' dx 


(continued) 
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9z l 

ay z 3 

<=fr> 

ay 

9z 2 

.2 

<-^f> 

dx 

ax z 4 

az 

ay~ = z 5 

( - a2w > 
' dxdy 

^3_ 

2 

(_ 8w ) 
' 3x3y ; 

ax z 5 

az 3 

ay z e 

<-T» 

ay 

!!i_ 

ax z 7 

_3 

(= 

1 3 ' 

ax 

az 4 

ay z s 

, 

' 2 ’ 
ax ay 

9z s_ 
ax z 8 

.3 

(= aw > 
dx xy 

9Z 5 

ay z 9 

a 3 w 

' 2 
ax ay 

az e 

ax z 9 

.3 

(_ 8w ) 
ax ay 

ay Z 10 

a, 3 ’ 

az 7 

ax u i 

( ,°s 

( ax 4 
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ay ' u 2 


(= 


4 

8 w 

ax ay 



u„ 


!!i _ 

ay u 3 


ax = u 3 


(= 


(= 


(= 


4 

aw 


a 

-) 

8x 8y 

„4 


a w 


2 

2 

ax ay 

ft 4 


a w 



2 2 ' 
dK By 


By u 4 


Bz 

10 

9x U 4 


8z io 

ay u 5 


(= 


(= 


(= 


4 

aw 


dxdy 

4 

8 w 
3x 3y 3 


4 

aw 
__ 4 

ay 


> 


) 



(=-) 
v dx' 



(=— ) 
1 3y ; 



k ^ 2 ; 
dx 


at 


— = [|3 (t x +K)2 1 -t 1 (u 1 +2» 3 +u 5 )-2t 2 (» 7 +V- 2t 3 (Z 8 +Z 10 ) " , l ( V VZ 6 ) 

(continued) 
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~V Z 6 +VZ 4 )]/2(1 ~ V)Z 5 


v 8x9y ' 


9x 3y 


{ 3x3y ' 


ay v 3 


<_ a/ ' 


(2.1) 


The boundary conditions, in the case of a rectangular boundary and simply- 
supported edges read: 

Z 1 
Z 1 


= 0. z 4 = 0 

along 

x = 0, a 


O 

II 

CO 

N 

o' 

11 

along 

y = 0 , b 

(2. la) 


We now form the Hamiltonian: 


H ‘ ‘i + W V-rtV W Va +x 6V : \V W S>VhoVhi*2 

+ Vl + WV5 + WV8 + V9 + VlO + WWW |i 10 U 5 


+tl ll t 3 ^13 V 3 2(l-v)z 


- 2t 3 (Z 8 +Z 10 ) - V l ( V VZ 6 )_V 3 (Z 6 + «4> 


f 2 

< (3 (t 1 +/c)z 1 -t 1 (u 1 


+ 2u 3 +U 5 ) “ 2t 2 (Z 7 +Z 9 ) 


( 2 . 2 ) 


The necessary conditions read: 


d \ _ 8H_ _ \d + ^12 2 

8x 6y “ 9z " 2(1- v)z ^ ( t i +K ) 

1 O 


^2 „_aH_ 

9x + 3y 3 z 2 ~ 1 


(continued) 
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fUx a + p =0 

Bug 8 ^7 


9H_ t 1^13 +H 'l2^ 

9u g " 9 + ^8 (l-v)z g 


= 0 


™-,X + p =0 
9u 10 ^9 
4 


t (A. +p ) 
9H = _ l v 13 12 ; 

9u c ~ ^10 2(1- v)z 

5 5 


= 0 


X + p 

9H , 13 *12 , , „ 

X (z + vz ) = 0 

9v 12 2(1- v)z 1 4 6 ; 

X i) 


X + LL 

9H _ 13 *12 

9v _ ^13 2(1- v)z 

u O 


( V vz 4 ) = 0 


(2.3) 


For a rectangular simply-supported plate, the boundary conditions read: 



along x = 0, a 



along y = 0, b (2. 3a) 

We therefore have to solve a system of 46 equations in the 46 unknowns 
z,u, t, v,X and p , which is a formidable task if we consider that the majority of 

r+-/ 

the equations composing it are partial differential equations. A logical and simple 
assumption, in the line of the one made for the shear plate problem (see 
Appendix l),will however yield a simple partial differential equation to be satisfied 
by the optimal displacement mode w. We verified that this assumption was indeed 
valid in the shear plate case, and similar ones also have proven to be valid in 
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14 

one-dimensional cases. We will come back to this assumption in length in 
Part E. 

The top equation of (2. 3) is similar in form to the original constraint 
equation (1. 10) which might be written under the form: 

8P 9Q 2 . 

— + = 6 (T +IC)w 

9x 9y ^ ' 


with 






This leads us to assume 




(Z 7 +Z 9 )+t 2 (z 4 + vz 6 )+(1_ 


(Z 8 +Z 10 )+(1 ' , ' )t 2 Z 5 +t 3 <Z 6 + 


v)t 3 Z 5j 
VZ 4>} 


7k. + n, Z, 

13 H-2 = + _J. 

2(l-v)z a 

t) 


(2.5) 


where a is a proportionality constant. Under this assumption, compatible with 
the boundary conditions when z x is prescribed to vanish along the edge, 


Z 1<V VZ 6> 

*12 = + 3 


Z 1 ( V VZ 4 ) 

1*13 ' + 3 


( 2 . 6 ) 


The three equations above (the control equations) are rewritten as 
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This relation simplifies into: 


2 2 2 2 2 
z. + 2vz.z +z +2(1- v)z == 0 ! + p z 
4 4 b b 5 1 


or, in terms of the sole variable w: 


' dx ' dx Z 9y 2 '8y 2 ' '•9x3y-' 


2 2 
w 


that we might rewrite as: 


2 2 P ^ 8 2 w \ 2 9 2 w 8 2 w1 ' 

(V w) + 2(l-v ) I I - — H=a;+(3 

»-V 8x8y' 8x 8y J 


2 2 
w 


( 2 . 10 ) 


The expression on the left hand side is, up to a constant multiplicative 
factor, the potential energy of deformation density of the sandwich plate, which 
is identical to the one given in Ref. 33 for an elastic plate with constant thickness. 
This potential density might also be expressed very simply as a symmetric 
quadratic form of the principal radii of curvature p and p of the deformed 
plate, and more precisely in terms of the so-called mean and total curvatures 
(ref. 33, p. 250). For a sandwich plate with varying thickness such as the one 
we consider for which the bending rigidity is simply proportional to the thickness 
t, the potential energy density is independent of t. This is actually not the case 
for a classical elastic plate with varying thickness h: as we will see in the next 
chapter, where the potential energy density will be encountered again, it is 
then a function of h, and more precisely is proportional to the square of the 
varying thickness. 

Actually, (2. 10) states for the sandwich plate an extremely general 
optimality principle encountered in a very broad class of structural optimization 
problems: for a given fundamental frequency of vibration, (2. 10) expresses the 
conservation of the difference between the potential energy and kinetic energy 
densities. We will come back in Part E to this extremely important feature, 
which is a necessary and sufficient condition for an extremal. 
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If we rewrite the left hand side as: 



we see that it cannot be negative, the Poisson ratio v of the given material 
having to be comprised between the extreme values 0 and 1/2. On the boundary 
3D, w = 0, and the above expression takes, from (2.10), the value a: this shows 
that a is actually a positive quantity; let: 


2 


a = c 


If w is a solution of: 


2 2 

(V w) + 2(l-v ) 



a/w aV) 

£K Z dy 2 J 


2 2 
= 1 +P w 


( 2 . 11 ) 


then ecw is a solution of (2. 10). As the displacement mode is known up to a 

multiplicative constant factor only, and if we recall that a was any multiplicative 

constant, it is sufficient for our purpose to consider (2. 11). The latter is a very 

interesting looking nonlinear partial differential equation of the second order in 

the dependent variable w, much simpler than what we might have expected from 

an original system of 46 equations. It is also very general in that it has to be 

2 

satisfied inside a domain D of any shape, the constant (3 being determined 
from the uniform plate, according to (1. 12). For a rectangular domain, 



The different boundary conditions corresponding to different supports (clamped 
or simply-supported edges) are to be chosen among (1. 10a). 

Once the optimal displacement mode is found from (2. 11), the optimal 
thickness distribution is given by the original constraint equation: 


96 



3. 


v 2 (r v 2 w)+(i 


-v)( 2 ^ 

' flva 


9 2 w _ 9 2 t a 2 w 

2 2 
9x9y 9x9y 9x 9y 


9y 9x 


^-yj-p 2 (r+K)w = 0 


( 1 . 10 ) 


A Numerical Application 

Equation (2.11) which gives the optimal displacement mode is, as 

already mentioned, a second order nonlinear partial differential equation. 

It is unfortunate but true that equations of this type have generated 

practically no interest among mathematicians in our century. The main 

39 

reason, as quoted from Ames seems that nearly no physical situations 
within the scope of mathematical formulation yields equations of the afore- 
mentioned type to be solved as steps towards the solution. It is hoped that 
optimization problems making use of variational techniques will give a new 
life to the whole subject: the constraint equation which we started from 
was a fourth order one, linear in w, whereas the optimality condition 
(2. 11) is a second order nonlinear one. Similarly, for the shear plate, 
the constraint was a second order partial differential equation linear in 
w, yielding an optimality condition in the form of a first order nonlinear 
equation. The characteristic features common to such problems seem 
that, whereas the constraint is linear, because based on a previous 
linearization of the system under investigation, the optimality criterion 
yields a nonlinear equation of order reduced by half. For sandwich struc- 
tures, we are always led to an equation involving the optimal displacement 
mode only'*''*’: it will be nonlinear of the second order if the original con- 
straint is of the fourth order, as is extremely common in the theory of 
elasticity. 

This present disinterest in nonlinear second order partial dif- 
ferential equations is the more curious because extremely interesting studies 
of the subject have been made at the end of the XYIIIth and beginning of 
the XIXth century by some famous mathematicians. In those earlier days, 
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the goal was then to find a general solution to the partial differential equation, 
and to impose the boundary conditions to determine the unknown functions or 
constants appearing in the expression of the general solution. This method of 
course does not apply to every case, but is worth of more efforts on the part of 
researchers confronted with a first or second order nonlinear partial differential 
equation. It worked in the case of equation (B. 2. 8) for the most general shape 
of the boundary. 

The author wishes to express his gratitude to Ames for referencing in a 
39 

modern work the admirable work by Forsyth, first published in 1890 and 

40 

reprinted by Dover in 1955 , which is a monumental tribute to the theory of 

differential equations and bears most of the knowledge in this field up to the 

beginning of the XXth century. In there were we able to find a mention of the 

original work of Ampere and in particular the exact reference of his two funda- 
41 42 

mental contributions ’ which led us to the discovery of an exact solution of 

equation (2. 11) when the boundary of the plate is an ellipse. However, we have 

been unable yet to find an exact solution of this equation for a rectangular 

boundary: the solution cannot in this case be simply a function of the distance to 

the boundary as it was for equation (B. 2. 8), due to requirements of continuity 

for the bending slope 0 , the bending moment and the shear along the 

lines loci of the points from which two different normals with equal length can be 

drawn to the contour, requirements which cannot be met under such an assumption. 

A numerical solution at this stage of research seems the most profitable 

thing to do, as it also will help orienting future investigations of equation (2. 11), 

The reference sandwich plate under consideration will be constructed as 

follows:* the face material is constituted of 7075-T6 aluminum alloy, which has a 

6 2 

Young' s modulus of 10. 4x 10 psi, a Poisson ratio of 0. 33 (1- v = 0. 89) and a 

3 

density of 0. 1007 lb/in . Its thickness is t^ = 0. 012 in. The core is made 

*The author wishes to express his gratitude to Professor Jean Mayers and 
Professor Richard Shevell for their help in getting the above data. 
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of 3/16-5052 H39-0015 P hexagonal aluminum honeycomb alloy, with a density 
3 

n of 0. 00254 lb/ in , and its thickness is h = 0.3 in. 
r c 0 

For a simply-supported square plate with sides of length a = 10 inches, 
the fundamental frequency of vibration is given by 


2„ 2 / 


(0.3) 2 x 10.4 X 10 6 


0. 012 


2x0. 89 


0. 3x0. 00254+0. 012x 0. 1007 


= 353 cps 


The values of the parameter k is 


■■ 2 s S is0ii! 


and p is given by: 


2 4tt 1 239 

(3 = - — — — — — ■ — = 0. 0239 
^ 1+K 4 4 

a a 


The numerical problem is now reduced to solving simultaneously the partial 
differential equations 


(V 2 w) +1.34(w 2 -w w ) = 1 + 0. 0239w 2 
xy xx yy 


(3.1) 


2 2 

V (r V w)+0.67(2t w -t w -t w )-0. 023 9(r + 0. 632 )w = 0 
xy xy xx yy yy xx 


(3.2) 


in the square domain limited by the straight lines x = 0,10 and y = 0, 10 under 
the boundary conditions (simply-supported edges): 


w = 0, 

T w = 0 

XX 

along 

x = 0,10 

(3.1a) 

w = 0, 

T w = 0 

yy 

along 

y = 0,10 

(3.2a) 
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For the computation, 49 x 49 = 2401 interior mesh points were considered, 
therefore breaking the domain into 2500 small squares with side length of 0.2. 

The initial start for w was taken to be either w = 0 or w = sin — sin the 
displacement mode for the original uniform plate. The solution was found from 
(3. 1) after the same number of iterations, 10, using a regular relaxation method. 
The contour lines for the optical displacement mode w are drawn in Fig. 14 
(interval taken equal to 0. 2). The maximum of w is attained at the center of 
the plate, and its value is equal to 1. 23. This pattern is also shown by the 
computer printout reproduced in Fig. 15. The letters A,B,C, . . . correspond 
to points with a w in the respective ranges 0. 1 ± 0. 015, 0. 2 ± 0. 015 
0.3 i 0. 015, . . . 

The above values for w are then used in (3. 2) which is now a second 
order linear partial differential equation for the optimal nondimens ional thickness 
r . 

Two initial guesses for r were taken, respectively, to be 

T= 1.2 sin — sin— and r =1.5 sin — sin 

10 10 1 5 5 1 

Both vanish along the edges, as should do the optimal thickness, because neither 
w^ nor w are zero along those edges. The convergence was in that case very 
slow, and the solution was found after 150 iterations. At that stage, the sum of the 
squares of the differences between the values of r during two successive iterations, 
sum taken over the 2401 mesh points, was only 0. 08736, thus securing the validity 
of the solution. Contour lines are represented in Fig. 16, while a computer 
printout of the optimal configuration, under those same rules governing the 
representation of w, is reproduced in fig. 17. The optimal thickness distribu- 
tion along a median axis of symmetry and along a diagonal of the plate are 
plotted infigs. 18andl9, respectively. The mass of the optimal face layers is 
69. 8% that of the original uniform ones, and this corresponds to a total mass 
saving of 18. 5 % for the sandwich plate. 

A computer program written in the FORTRAN H language to solve the 
system (3. 1), (3.2), using double precision, is presented in Appendix 3. 
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r vanishes only along the edges, which does not prove to be practically 
inconvenient, as was the optimal shear plate. Also, its maximum is 1.47, which 
is quite reasonable. However, one might be interested in imposing a minimum 
thickness constraint as an extra inequality constraint. The application of the 
general method outlined in Part A is straightforward, and the above procedure 
may be used with only slight modifications. Due to limitations in computer time, 
numerical computation was not carried out, but it is a simple exercise for the 
mind to imagine how the optimal thickness distribution will look like when an 
inequality constraint is applied. 
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Figure 14. Contour lines for the optimal displacement mode w-square 

2 

simply- supported sandwich plate (the constant c has been 
set equal to 1). 
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Figure 15. Computer printout for the optimal displacement mode. 

Letters A, B, C, . . . correspond to points in the range 
0. 1 ± 0. 015, 0. 2 ± 0. 015, 0. 3 ± 0. 015, 
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Figure 16. Contour lines for the optimal thickness distribution 
square simply- supported sandwich plate. 
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Figure 17. Computer printout for the optimal thickness distribution. 
Same notation as in fig. 15. 
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1.2 



Figure 18. The optimal thickness distribution along a cross-section 
following a median axis of symmetry. 
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PART D 


THE MINIMUM-MASS DESIGN OF AN ELASTIC PLATE FOR A GIVEN 
FUNDAMENTAL FREQUENCY OF VIBRATION 


Derivation of the necessary conditions leads for vanishing displacements 
along the boundary to an optimality criterion expressing again the conservation of 
the Lagrangian density. It differs from the conditions previously encountered 
for sandwich structures in the fact that the design parameter h, thickness of 
the plate, is present. The optimal displacement mode and the optimal thickness 
distribution are therefore found to be the solutions of a system of two simultaneous 
nonlinear partial differential equations. Only a numerical approach seems 
possible in that case. 


1. Statement of the Optimization Problem 

The differential equation of equilibrium for a plate of variable thickness 

43 

h can be found in the classical textbook by Timoshenko, pp. 173-174 or in a 
44 

paper by Reissner . In rectangular coordinates (x,y) in the plane of the plate, 
it reads: 


V 2 (DV 2 w)+(1-v) 


Uis. 

' 9x9y 


9 2 w 

9x9y 


9 2 D 9 2 w 
9x 2 9y 2 


9 2 D 9 2 w \ 
9y 2 9x 2 / 


in which 

w denotes the lateral deflection 


D = — the flexural rigidity 

12(1- v) 

v the Poisson's ratio and, 
q the intensity of the acting load. 


( 1 . 1 ) 
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The free vibration of such a plate is governed by equation (1. 1) in which 
q is replaced by the d' Alembert force 


q = -ph 



( 1 . 2 ) 


p being the density of the material, assumed to be homogeneous, of which the 
plate is made: p is therefore a constant. 

The free vibrations of the uniform reference plate made of the same 
material and of constant thickness h Q are governed by the equation 


n _4 9 2 w 

D V w--ph ~J 

at 


where 


D o = 


Eh o 

12(1- v 2 ) 


For a simply-supported rectangular plate extending over a domain D 
defined by 0 ^ x ^ a and Osysb, the natural frequencies of the system are 
found to be 


GJ 

mn 


2 

TT 




Eh, 


12p(l-v ) 




(1.3) 


The corresponding natural modes are 


mTrx mry 

W (x,y) = A sin sm — rf- 

mn mn a b 


(1.4) 


The fundamental frequency of vibration is 

O ' 


W f °11 


” 2 (Vj) 

' Q h » 


Eh 


0 


(1.5) 


a b 'V 12 p( 1 ~ v ) 
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The variable thickness h of a plate occupying a domain D and having 
the same fundamental frequency of vibration co^ than the uniform one made of 
the same material and occupying the same plane domain D under the same 
boundary conditions, together with the corresponding mode w, have to satisfy 
the partial differential equation: 


V 2 (DV 2 w)+(l-v) 


( 2 ^ 
' avflv 


2 2 2 2 2 2 
" aw 9 D 9 w 3D " 


dxdy dxdy 3x^ 9y 2 8y 2 dx 


aM 


-phw^w = 0 


( 1 . 6 ) 


For a simply-supported rectangular plate, the value of is given by 


(1.5). 


The constraint equation (1.6) may be rewritten in terms of h and w 


only, as 


24 [ah a 2 ah a 2 

h V w + 6 h[^— -(V w) + --(V w) 


+ 3 


3 


C [*(S? •*$]&••?) -‘"•"[■i !»&]& 

Gfe)' •■&]&*•&)} 


9y ' 9y dx 

or, as an alternative way, in terms of D and w only as 


12(l-v ) 2 

— ^ P w f w =° (1-7) 


( z z z 

a d a w 3D 

2 ~ - — £ 

rbcrtv dxdv fW 


2 d a 2 w a 2 D a 2 D^ 21/3 


8x3y dxdy dx“ 9y 2 9y 2 9x 2 


)- 


P D w = 0 

( 1 . 8 ) 


where the positive constant p is defined as: 

1/3 


P 2 = 


[^J 


P 


(1.9) 


Our goal is to minimize the total mass of the plate, given by the surface 
inter gra.1: 



M = J'J’ p hdxdy 


D 


and the optimization problem may be stated as follows: 
Minimize the functional 


( 1 . 10 ) 



D 


subject to the constraint (1.7), or, alternatively, 
Minimize the functional 


( 1 . 11 ) 


K 


ff D l/3 dxdy 


( 1 . 12 ) 


subject to the constraint (1.8). 

The boundary conditions can be expressed for the more general domain 

D in terms of the coordinates n and s in the directions normal and tangential 

32 

to the boundary respectively: 

At a simply-supported edge, 

w = 0 

2 / 1 9w 9 2 w \ 

-DV w ♦ (l-v)D^ - _ + — j - 0 


For a clamped edge, 
w = 0 



9n 
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In the case of a free edge, 


-DV w +(l-v)D 


(l =0 

VE8n q 2 ) 


n 9 ,2 , 9 r_/ 9 2 w 1 ten A 

D an (v W)+( ~ V 9s L \ 9n9s “ B 9s j J 


where the Laplaeian has the form 


v 2 = — i 1 9 i Q 2 

9n 2 R 8n 8s 


and R denotes the radius of curvature of the boundary curve 3D. 

For a simply -supported rectangular plate, the boundary conditions read: 


w = 0 

and 

a 2 
o w 

D-~ 0 
ax 2 

along 

x = 0, a 


w = 0 

and 

fl 2 
O w 

D-f = 0 
9y 2 

along 

y = 0,b 

(1.8a) 

For clamped 

edges, they read: 




w = 0 

and 

^ = o 
ax 

along 

x = 0, a 


w = 0 

and 

9w 
9y " 

along 

y = 0,b 

(1.8b) 


2. The Necessary Conditions 


We will concentrate on the optimization problem under the second 
formulation, where we wish to minimize the surface integral 


K 


ff » 1/8 d*dy 


( 2 . 1 ) 
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subject to the constraint (1. 8), rewritten as: 


_ _4 3D 9 2 „ n 3D 9 ,2 . 

Dv w + 2 ST^ (V w) + 2 iFiF (V w) 


d 2 p f 9 2 w 9jv\ . 9 2 D 8 2 w 3 2 D / 9 2 w d%,\ 

2\ 2 + V 2 ) + (_V> 9x9y 9x9y 2\ 2 2J 

dx x dx dy 7 dy x 9y 9x 7 


- P V /3 w = 0 


( 2 . 2 ) 


In a similar fashion to that which was done in Part C, we break (2. 2) into 
a system of first order partial differential equations as follows: 


z^ = w 


9z, 


= z„ 


9x 2 
3z, 


9y 3 

ff2 

3x Z 4 


(=— ) 
' 3x ' 


(=— ) 
' 9y ; 


(= 


9 2 w , 

9x 2 


!!2 

9y z 5 

!fs_ 

9x Z 5 

!!3 

dy Z 6 


a 2 
9 w 

9x9y 


a 2 
9 w 

dxdy 


2 

3 w. 


dy 


9x Z 7 


(=^) 

' „ 3 ' 

3x 


(continued) 
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<= — ) 
v 2 
dx dy 


,3 

(=4“) 

9x 9y 


a 3 

(= ^_) 
k 2 
dxdy 


a 3 

(=5-«L_) 
1 2 
dxdy 


J 

(=— ) 
' „ 3' 
9y 


<- 4 > 

9x 


9y U 2 


(=4-) 

dx°dy 


9x U 2 


4 

(= § -| L - ) 
9x 9y 


9y U 3 


4 

' „ 2 „ 2 ' 

9x 9y 


9x U 3 


4 

v 2 2 ' 
dx dy 


dy U 4 


<=^T> 

dxdy 


(continued) 



3z 


10 

9x - U 4 
3z. 


10 


u_ 


3y 5 


^=D 
3x 2 


(= 


4 

3 w 
3x3y" 


3 4 w 

<“— > 

9y 


^=D 

9y 3 


^2_ 

3x V 1 


(=— ) 
* n 2> 
3x 


9D 2 2 1/3 

— = [P D z 1 -D(u 1+ 2u 3 + u 5 )-2D 2 (z 7+ z 9 )-2D 3 (z 8+ z 10 ) 


- y 1 (z 4 + v z 6 )~ v 3 ( z 6 + vz 4 )]/ 2 ( 1 " v ) 1 z 5 

3D„ 3D„ 


3x 3y 


v 3x3y ’ 


9D 3 

9y v 3 


(= a2 °) 
v 2 ’ 
9y 

(2.3) 

In the case of a. rectangular boundary with simply-supported edges, the 
boundary conditions read 

z_ = 0, z = 0 

1 4 

along 

x = 0, a 


Z l = °’ Z 6 = ° 

along 

y = 0,b 

(2.3a) 

We form the Hamiltonian: 
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H 4 D 1/3 4 ^ + V 5 + V 7 + Vs + V« + Vl + Vz + W Vi 


+ \l D 2 + \ 2 V “lV ' i 2 Z 5 + V 6 + “4 V V 9 + VlO + »7 V V 3 + Vi 
+ ‘ i 10V“ll D 3 +l *13 V 3 


A 13 +Ml 12 

2(l-v)z 


^P 2 D 1/3 Z 1 -D(n 1 
- v 1< Z 4 +VZ 6 ) - V 3 <Z 6 +VZ 4 ) } 


+ 2u 3 +U 5 ) - 2D 2 (Z 7 +Z 9 ) - 2D 3 (Z 8 + V 


The necessary conditions are: 


9A. 0m, X + n. . 

1 + __1__9IL 13 ^12 2 1/3 

3x + 9y 9z^ 2(1 -v)z 5 ” 


^2 + ^2 

3x 9y 


|2-=-X 
9Z 2 1 


az 3 -i 


9H A 13 +U 12 

3z 4 " A 2 + 2(l-v)z g (v i +vw 3 ) 


9H , A 13 + ^12 

= — A, —M, H 

8z* 3 2 0/1 X 2 

5 2(l-v)z 

5 


|^P 2 D 1/3 z 1 -D(u i 


+ 2 V U 5> 


- 2D 2 (Z 7 +Z 9 ) - 2D 3 (Z 8 +Z 10 ) - V 1 ( V VZ 6 ) - V 3 (Z 6 +VZ 4 ) 


9H 13 *12 , 

3z ~ _U 3 + 2(1- v)z 5 (V 3 +VV 1 } 


9H , A 13 + ^12 „ 

— = -X + D 

9z 4 (l-v)z 2 


(continued) 
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For a rectangular simply-supported plate, the boundary conditions read: 


z = A =A =z =A =A = A =A =A =A = A =o 

1 2 3 4 5 6 7 8 9 10 11 12 13 


along x = 0,a 


2 l“ ‘V V V V V V V V 'V >*11" ‘‘l2 = ‘‘is = 0 


along y = 0,b 


(2. 5a) 


As in the sandwich plate problem of Part C, we note that the first 
necessary condition 

3A 8ij. A +u. _ „ . 

1 _1 = _ 13 *12 2 1/3 

ax ay 2(1- v)z P 

O 


is very similar in appearance to equation (2. 2), rewritten as: 

|£ + = pV'V 

3x 3y 1 

where 



and therefore leads us to assume the proportionality relations: 
V ! = [D(2 7 +z 9 ) + D 2 (2 4 + vz 6 ) + (1 -v)D 3 z 5 J 

Q 1 

u, = — = - - [D(z +z, )+(l-v)D„z_ +D„(z +vz.)] 

P 1 a a 8 10' ' ' 2 5 3 V 6 4 /J 

A +n, z 

13 *12 _ _1 

2(l-v)Zg a 


(2.6) 


(2.7) 
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or, in terms of the original w and D variables: 


2 2 f q2w ^ 9 2 w 9 2 w1 

(V w) +2(1- v) I I 1 — 

l-'dxdy' 9x dy J 


1 ~-2/3, 2 2 

= -D (a + (3 w ) (2. 10) 


Equation (2. 10) is a second order nonlinear partial differential equation 
in the two dependent variables D, optimal flexural rigidity, and w, optimal 
displacement mode. Together with equation (1. 8) it forms a system of 
simultaneous partial differential equations leading to the solution of the optimiza- 
tion problem. 

(2. 10) can also be rewritten in terms of the optimal thickness distribution 
h and w, as: 


2 2 f7 8 2 w 'f 9 2 w a 2 w"l 

(V w) + 2(1- v) I ( -^-1 l = 

LAaxay/ ^ q/J 


2 2 
4(l-v ) 2 a!*+w 

“ E P “f ~2~ 


( 2 . 11 ) 


the constant or* being defined as 


O’* = 


a 



and, together with (1. 7), forms a system of two simultaneous partial differential 
equations for the two unknowns h, optimal thickness distribution, and w, corre- 
sponding displacement mode, the latter being known up to a constant multiplying 
factor only, as usual, a* can be shown to be positive, exactly as for the sand- 
wich plate equation. Let 


4(1-^) 

E 



a* = 


2 

c 


The solution w will be of the form cw^, where w^ satisfies: 

2 2 r f 9 2 w 'f 9 2 w 9 2 w"1 4(1- v 2 ) 2 2,2 

(V w) +2(l-v) ( ] - — — — = 1 + -^-r — 1 pu w /h (2.12) 

L ^9x8y/ 9x 2 8y 2 J E f 
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For a simply-supported rectangular plate, the solution of the optimization 
problem reduces to that of the two simultaneous partial differential equations: 


K 2 \2 2 2 q 4 9 

J ““TT = ~ 

9x9y ' dx 9y J Sab 1 


ce* +w 


(2. 13) 


2 4 

r V w + 6 T 


[~9t a_ 

|_9x 9x 


2 9r 9 2 

(V w)+ i7 w) 


0 


a- 2 *2 

+3 l 2 <s ,+T ri 

9x 


« 2 A 2 

d w d w . 

+ V — 1 + 2(1+ V) 


L 9x 


9y 


r„ 9 t 9r + 9 2 T 9 2 w ~ ] 

9x 9y 9x9j 


9 r 2 9 2 T 

2 ^ )+r ri 


9y 


9 w 9 w ] 

— + v — — 
2 2 
k 9y 9x 


4 1 1 2 

w = 0 

a b 


9x9y 9x9y J 


(2. 14) 


T = h/h^ being the non-dimensional optimal thickness distribution (r = 1 for the 
uniform reference structure). 

The boundary conditions for simply-supported edges are: 


w = 0 

and 

9 w . 

T A 2 

dx 

along 

x = 0, a 


o 

II 

is 

and 

9 2 w_ 

T -J-0 

9y 

along 

y = 0,b 

(2. 14a) 

For clamped 

edges: 




o 

II 

£ 

and 

f£ = o 

dx 

along 

x = 0, a 


© 

II 

£ 

and 

9y 

along 

y = 0,b 

(2. 14b) 


A computational approach at this stage of the optimization process seems 
now inevitable, and the above system seems suitable to such an approach. 
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Equation (2. 12) is very similar to the optimality condition derived by 
45 

Suhubi for the optimal plastic design of a plate using an entirely different 
approach. 
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PART E 


THE OPTIMALITY CRITERION: A SUFFICIENT CONDITION 
FOR AN EXTREMAL 


The nonlinear partial differential equation to which the entire system of 
necessary conditions reduces in each of the three cases considered so far under 
the assumption of prescribed vanishing displacement at the edge is the expression, 
in each case, of an optimality criterion expressing the Constance of the energy 
distribution of the system, which is a sufficient condition for an extremum and 
was first derived by Prager for a broad class of structural optimization problems. 
The system of two partial differential equations to be satisfied simultaneously or 
not by the optimal displacement mode and the optimal thickness distribution, 
which is formed by this optimality criterion and the original constraint equation, 
and to which each of the aforementioned cases reduces is therefore necessary and 
sufficient. For the shear plate, uniqueness of the optimal solution follows. 

When the displacement is not required to vanish along the edge, Prager' s 
criterion still applies, but has to be written in a three-dimensional space. The 
above fact agrees with the validity of the proportionality assumptions made in the 
course of the derivation of the optimality condition. 

* * 

Equations (B. 2. 8), (C. 2. 10) and (D. 2. 11) are nothing but the expression 
of a very general optimality criterion stating the uniformity of the energy 
distribution throughout the structure under consideration, valid for three dimen- 
sional structures to be optimally designed under the most general assumptions. 

46 47 

It was first derived by Prager ’ and is an alternate form of stating the so- 
called "conservation of the Lagrangian density" which was pointed out by Ashley 
and McIntosh.^ 


t 13 

The term Lagrangian density which was introduced by Ashley and McIntosh in 

optimal frequency design has its origin in the similarity between the Lagrangian 

of Classical Mechanics, difference between the kinetic and potential energy, and 

the expression G-cp H, difference between the potential energy density and the 

quantity pw^w^, which can be viewed as a kinetic energy density. 
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46 

Using Prager' s notation, the structure is to be designed to satisfy 
a single behavioral constraint prescribing the value of a scalar $ itself 
characterized by a global minimum principle involving a field cp associated with 
the structure, and having one of the following forms: 


$ 


$ 




min 1 F[cp]dV 


V 




G[cp]dV 


mm 


V_ 

/ 

V 


H[cp]dV 


(la) 


(lb) 


V being the volume occupied by the structure. When the constraint is on the 
fundamental frequency of vibration, is the square of the prescribed value 
Uj of the frequency and, from Rayleigh' s principle follows that the constraint is 
of the form (lb) above, where G[cp] is twice the strain energy density of the 
field cp, and H[cp] is p[cp| , where p is the density of the material. 

Prager has proven that when the behavioral constraint is expressed by a 
minimum principle of the form (la), a sufficient condition for an optimum is: 

F = Constant (2 a) 

on the portion of the surface S where vanishing surface tractions are prescribed. 
When the constraint is in the form (lb), a sufficient condition is: 

G - $ H = Constant (2b) 

on that same portion of the external surface S. 

The above conditions can also be shown to be necessary. * 

*The general formulation of the sufficient conditions above is actually slightly 
more complex, and the best way to express them is, as did Prager, to use 
notations from set theory. The reader is referred to references 46, 47 for more 
precision. 
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A proof of the sufficiency of the condition of uniform energy distribution, 
which we believe to be original, will be given below for the case of interest to us, 
the minimum-mass design of a structure (one, two or three-dimensional) for a 
given fundamental frequency of free vibration. This proof makes use of Rayleigh' s 
principle. 

Let & be the set of all structures with the same constant density p 
similar to the uniform reference structure D Q and all having the same funda- 
mental frequency of free vibration co^ under the same support conditions. 

Among those, at least one, which we shall call D, is assumed to be such that its 
energy distribution is constant throughout the volume V it encloses. If cp is 
the displacement field associated with it, G[cp] twice the strain energy density 
of the field, we assume that: 


2 2 

G[cp] - w f pep 


2 

c 


(3) 


and, from Rayleigh' s principle, is exactly given as: 


f G[ cp ]dV 


2 V 
“f 


J' pcp 2 dV 


V 


rewritten as: 


(4) 



(4a) 


Now consider any structure D member of the above set and occupying a volume 
V, with the corresponding field cp. From Rayleigh' s principle, 
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( 5 ) 


J G[cp]dV 

2 V 

w f = 

/pcp 2 dV 

V 

and, cp being the field corresponding to the structure D, we also have the 
inequality: 

(6) 

V 



rewritten as: 



(6a) 


Subtracting (4a) from (6a), 


J SG[cp] - c^pep 2 J 

v-v ^ ^ 


dV^O 


and, by the use of (3), the inequality reduces to: 


( 7 ) 


/ 

V-V 


dVs 0 


or 

V <; V (7a) 

The volume V occupied by the structure D satisfying the condition (3) 
is therefore smaller or equal to the volume V occupied by D. As the structure 
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D was arbitrarily chosen among the set $, we are led to the conclusion that a 
structure satisfying (3) has the smallest possible weight*' among the structures 
belonging to the class & under consideration. (3) is therefore a sufficient 
condition for an extremal of the weight under the conditions already stated. 

The above sufficient condition is intuitively not at all surprising: this 
uniformity of the energy distribution is a very general and logical property 
encountered whenever one wants to make the best out of a given system, and 
should be considered as a law of nature. In the case of a vibrating plate with 
uniform thickness, the efforts are not uniformly distributed throughout the plate, 
and some parts undergo more stress than others. We might say that the uniform 
plate is not ideally designed for vibration purposes, and it reacts to vibrating 
conditions in an unorderly and anarchistic behavior, which is however the best it 
can do if it cannot modify its total shape. This is also a law of nature. Suppose 
now that the plate could modify its shape and mass: it would still comply to the 
laws of nature, but with more freedom now, and tend towards an optimal shape. 
The non-uniform plate would now be perfectly adapted to the vibrating conditions 
imposed on it, and the efforts would now be equally distributed throughout. This 
is the best of all plates under the given conditions, such that the energy distribu- 
tion is uniform and the same role is assigned to any point of the structure. 
Understandably, such a structure will be optimally designed for the purpose at 
hand, here the first vibration frequency being held constant, and this optimal 
plate will have the lesser possible weight of all structures having the same 
fundamental frequency, its matter being the more properly distributed of all of 
them. 

This uniformity of energy distribution attained by Nature in the best of 

all configurations of a system is also a very general principle going beyond 

48 

structural optimization: Max Munk showed in 1918 that the minimum induced 
drag of a wing is obtained if the distribution of lift over the span is elliptic, 
corresponding to a uniform energy distribution over the wing area. 

*The material is assumed to be isotropic with uniform density p. 
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Let us now verify that we indeed obtained the optimality condition for the 
three structures we considered. 

For the shear plate, the strain energy density is: 


v = 

2 u 8x 7 


aw .* 5 

+ W 1 


and H is simply: 

2 

H = pw 

so that (2a) is expressed under the form: 


( 8 ) 


( 9 ) 


,9w.^ ,9w 2 p 2 2 

<i? + V " o “f w " 00,lst,mt 


( 10 ) 


which is nothing but equation (B. 2. 8) which is therefore a necessary and sufficient 

condition for an extremum. This equation was derived for simply-supported 

edges, and was found to have a unique solution, given analytically by (B. 2. 9). It 

follows that it is the solution to the optimization problem, and the shear plate 

problem possesses a unique solution as described in Part B. 

The strain energy density of a sandwich plate is, in terms of the lateral 
46 

deflection w(x,y): 



expressing the optimality condition (2a ) as: 



2 \ 2 / 2 \ 2 
3 w \ „|9w 

+ 2 I 1 + 2v 

\3x3 yj 


2 2 
3 w 9 w 

2 2 
3x 9y 



2 2 

-(3 w = constant 
(12) 


128 



where (3 is a constant whose value is the same as the (3 defined by (C. 1. 12). 

This is nothing but another way of writing (C.2. 10). Our proportionality 
assumption (C. 2. 5) was thus perfectly valid, and (C. 2. 10) is not only a necessary, 
but also a sufficient condition for optimality. 

The fact that for sandwich structures the optimality condition does not 
involve the design parameter j was pointed out by Prager and Taylor. ^ 

For the classical plate with varying thickness, the strain energy density 
in bending is given by:^ 2 


V 


Eh 

8(l-v 2 ) 


2 2 
(V w) + 2(l-v) 





(13) 


and Prager' s optimality criterion leads to the expression: 


(V 2 w) 2 +2(l-v) 



n 2 .2 
9 w 8 w 

2 2 
9x 9y 


4(l-v 2 ) 22 

— ^ — 1 p w = constant 

(14) 


which is exactly (D. 2. 11). 

The precision and beauty of Prager' s theorem is even more evident when 
one considers the following important fact, which will also serve the purpose of 
clarifying a delicate point, which is the assumption of proportionality in the 
course of the derivation of the necessary conditions. This remark led to the 
conservation property; however, it was compatible with the boundary conditions 
only in the case of vanishing displacements at the edge: this was the case for the 
simply-supported shear plate, for the clamped or simply-supported sandwiGh 
plate and also for the elastic plate. On the other hand, Prager' s theorem 
seems to lead to the same optimality condition, whatever the boundary conditions 


129 



are. A careful review of the conditions under which it is valid is therefore 
necessary at this point. 

Prager' s theorem applies to optimal structures occupying the volume V 
with the surface S. Each element of S is supposed to belong to one and only 
one of the sets S' , S", S'", where nonvanishing surface tractions are prescribed 
on the elements of S' (for problems that we consider here S' =0), vanishing 
surface tractions on the elements of S" (here the surface covered by the plate), 
and vanishing displacements on the elements of S'" . Moreover, the surfaces 
S' and S'" and the boundary conditions thereon are regarded as fixed 
"ingredients" of the design problem. S'" is the boundary of the plate when the 
displacement is prescribed to vanish along the edge. However, the portions of 
the boundary which are not supported belong to S" , the plate being 
viewed in the broad sense of a three-dimensional structure, and cannot be 
regarded any more as fixed ingredients of the problem: in such cases, the 
optimization problem needs to be viewed in three dimensions , and will not yield 
an optimality condition as simple as (B. 2. 8), (C. 2. 10) or (D. 2. 11). 

Prager' s theorem therefore still applies, but with an added dimension. 

As we recall, when portions of the edge are free, the proportionality assumption 
expressed by (B. 2. 7), (C. 2. 5) or (D. 2. 7) is no longer compatible with the boundary 
conditions, rendering equations (B. 2. 8), (C. 2. 10) or (D. 2. 11) no longer valid. 

This is one more evidence of the complete agreement between Prager' s 
optimality criterion and our results, obtained by an entirely different method 
based on the calculus of variations in two dimensions. 
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CONCLUSIONS 


In our approach to optimization problems in two dimensions, outlined in 
the theory presented in the first part of this work and applied to the subsequent 
structural optimization problems, we are confronted with a system of partial 
differential equations to be solved inside a domain, together with boundary con- 
ditions, both expressing a set of necessary and sufficient conditions for an 

optimum. A computational approach undertaken at this stage, although being in 

49 

theory a classical problem of Numerical Analysis would prove very awkward 
and especially time-consuming: for a simple structure such as a sandwich plate 
or an elastic plate and under a very simple constraint such as the fundamental 
frequency of vibration being fixed, the necessary conditions yield a system of 
46 equations in 46 unknowns, mostly first-order partial differential equations! 
However, and this is one of the purposes of the present study, a careful study of 
these necessary conditions in each particular case permits, under assumptions 
later verified, to reduce them to only two partial differential equations in the 
important unknowns. In some cases we will be lucky enough to find an exact 
analytical solution, and in the majority we will have reduced the original problem 
to a reasonable amount of numerical computation, following the well-established 
pattern to solve numerically a single partial differential equation or a system of 
a limited number of the same. 

This property leads us to make the remark that, in the case of optimization 
problems in one dimension, the necessary conditions for an extremum lead to two- 
point boundary value problems, requiring extreme caution: no universal numerical 
method has been found up to date to solve any kind of such problems (for example, 
the powerful transition-matrix procedure which gave excellent results in a number 
of cases in structural optimization failed for the panel-flutter case where the non- 
linearity was too strong). However, for two-dimensional optimization procedures, 
following the method presented here, we are in the end confronted with a problem 
of a classical type and there is, in theory at least, a method permitting to solve it. 
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This is a very important feature differentiating optimization problems 

in one -dimensional space from those in higher dimensional spaces. This 

50 

is why we will be less pessimistic than R. Chattopadhyay about the 
existence of practical solutions for problems of optimal control of 
systems with distributed parameters, always keeping in mind that asking 
for the services of a computer should be the last thing to do, and should 
be done only after we have made sure that every possible method or 
trick failed. It is not our goal to discredit the computer: its help is 
tremendous in the majority of cases, and it brought the solution to 
inextricable practical problems, but some reflection makes it more ef- 
ficient and . . . cheaper! 

However, we will at this point raise several questions, as for- 
mulated in Ref. 50. 

Is it better to discretize the system itself in the first place, as 
19 

did Haug, than to discretize the equations expressing the necessary 
conditions for an extremal, obtained for the exact system? Are the ap- 
proximate equations consistent with the original system equations? Is 
the discretization of the system valid? Do the approximate solutions 
converge to the actual one in the limit and do the errors made remain 
bounded? These questions are unfortunately very difficult to answer 
and are beyond the scope of the present investigation, but on their solu- 
tion rests the possibility of obtaining results of practical interest. 

It is our sincere hope that the present study will lay the theoretical 
foundations necessary to future investigations in this entirely new field 
of two-dimensional structural optimization, reducing any such problem 
to the search of the numerical solution of a simple system of partial dif- 
ferential equations: this is by itself a classical but unfortunately non- 
trivial problem, and attaining the solution will rely mainly on the skill 
of the researcher as well as the special nature of the structure investigated. 
The choice of two-dimensional systems to be investigated is of course 
very wide, and more general constraints can be considered, still fol- 
lowing the pattern presented before. 
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APPENDIX 1 


INDEPENDENCE OF THE FORM OF THE SYSTEM OF PARTIAL 
DIFFERENTIAL EQUATIONS ON THE OPTIMIZATION PROBLEM 


In Part A, we showed that the constraint defining the optimal problem for 
the shear plate, expressed in the form of the single partial differential 
equation of order 2 in w and 1 in h: 



9w . 9 , 9w . 2 

— ) + — (h — ) + k hw = 0 
Bx ' By ' By 


(1) 


could be put under the form of two systems of very different appearance, part 
of the control variables being in one of them the thickness h, in the other one 
its first derivatives. 

Derivation of the necessary conditions for an optimum was made for the 
first case in section B. 2. 1 when the functional to be minimized is the integral 
of h over a plane domain D under the given constraint. Let us show now that 
the same answer is found if we start from the other system, obtained by the 
means described in (A. 1), and derive from it the necessary conditions. 

With this system as a constraint, the Hamiltonian takes the form: 


H = h + 


A z 
1 2 


+X 2 VW X 4W 3 


+ V2 


? 19 2 1 

^3<V k Z 1 + — > + ^2 < 2 > 


The necessary conditions for an extremal are thus found to be: 


d \ 

9 M 'i 

. 2 

(3) 

Bx 

+ 9y 

H 

oF 

* 2 

^2 

1 1 

- A Z 3 " 1 

(4) 

Bx 

9y 

" 1 + h 

CO 

% 

^3 V 2 

(5) 

Bx 

+ 9y 

-^1 + h 
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( 6 ) 


^4 1 V 1 Z 2 +V 2 Z 3 

9 x + 9y ^3 2 


V^3 = ° 


V^2 ~ ° 


(V 

(8) 


u„z„ 

x 

4 h 


(9) 


Vs 


^4 h 


= 0 


( 10 ) 


The form of the first equation suggests to try the substitution: 

hZ 9 

X. =• 


1 a 

hz 3 
X = — — 

2 or 


hz. 


or 


(ID 


by analogy with the original constraint equation. Now, 


x ._Vs 

4 or 


(9) 


V 3 

or 


( 10 ) 


and (6) is rewritten as: 


+ i<Va> ■“ - r WYs' 


( 6 ) 


or, in terms of w and h: 
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( 6 ) 


,9w 2 ,8 w 2 _2 w . 8h 8w 8h 8w. 

K dx > '8y ' h '3x 8x 9y8y ; 


But, from the original constraint, 


9h 8w + 9h 9w 
9x 9x 8y 3y 


2 2 
- k hw - hV w 


so that (6) is put under the form: 


,8w. 2 .9w. 2 .2 2 

W +( iF> ““ +kw 


( 6 ) 


which is nothing but equation (B. 2. 8), the solution of which led to the optimal 
thickness distribution of the shear plate. It is easy to check that the hypotheses 
(11) are consistent with the system formed by the original constraints, the necessary 
conditions and the boundary conditions in the simply-supported case. 


135 



APPENDIX 2 


ON ANOTHER METHOD TO SOLVE THE PARTIAL DIFFERENTIAL 
EQUATION ARISING IN THE SHEAR PLATE OPTIMAL PROBLEM 


A very interesting alternative approach to solve the equation: 


A 2 +{ SSL) 2 . 

K dx' K dy’ 


2,22 
c +k w 


( 1 ) 


and also applicable to more general boundary conditions than the one we considered, 
i. e. , w = 0 along 8D, has been suggested by Dr. E. O. A. Naumann. 

Let us make the change of dependent variable, introducing the new function 
z defined by: 



sinh(kz) 


(1) reduces to: 


9z 2 8z 2 

W + <ij> = 1 


or, using standard notations , 


( 2 ) 


( 3 ) 


2 

P 



( 4 ) 


This partial differential equation is a classical one, encountered in 
particular in the theory of plastic torsion of cylindrical beams . It also arises 
in geometrical optics, and is investigated under that name in Ref. 35. 

The system of five ordinary differential equations for the characteristics 
assumes the form: 



dy = 

ds 


2q 


(continued) 
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( 5 ) 


i = 2(p 2 + ? 2 > 


0 

ds 


fl = 0 

ds 


It follows that p and q are constants along any characteristic, and that: 
x - Xq = 2ps 
y - y n = 2ps 


z 



( 6 ) 


Hence the characteristics are straight lines. In view of the restriction (4) 
on p and q, they are precisely those lines which make an angle of 45° with the 
z-axis. A geometrical interpretation of those is given in Ref. 35, pp. 41-43. 
Elimination of s in (6) results in: 


P = 


x-x 


0 


z-z. 


q = 


y-y 

z-z 


0 

0 


and a solution of special interest is the cone of equation: 

(z-z Q ) 2 = (x-x 0 ) 2 + (y-y Q ) 2 (7) 

It can be used to form a complete integral of equation (1). To satisfy 
the boundary conditions, the following conditions need to be met: 

z = £> for x = % and y = q 

and dL = 0 along the given boundary curve t| = t|(§). 

Thus: 

(z— £) 2 = (x— £ ) 2 + (y-t) ) 2 (8) 

and, by differentiation: 
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(z-L)dC = (x-£)d£ + (y--n)dri = 0 


or 

(x-£) + (y-q)q' = °> T l' ( 9 ) 

Equations (8) and (9) together with the relation between | and q 
expressed as: 

q = q(£) (10) 

and its derivative with respect to | : 

q ' = V(£) (ii) 

are sufficient through elimination of § , q , q 1 to establish the desired function 
z = z(x,y) which satisfies the partial differential equation (3) and the boundary 
condition. For complex functions q(£), it may not be possible to find a single 
equation and % and/or q may remain as parameters of a set of equations. 
However, z is still uniquely defined. 

In the simple case where the boundary 9D is circular, the relation 
between § and q takes the form: 


1 2 2 2 
S +q =B 


together with the derivative: 

£ + qq' = 0 

Elimination of q ' from (9) and (11) gives: 
t|X= iy 

Substituting for £ into (8) and (10) gives: 


(z-f,) 2 = (x + y 1 ) + (y-q) 2 = (x 2 + y 2 )(l - ^) 


2 x 2 q 2 2 2 2 a 2 

E = + q = (x +y )(^) 


( 10 ) 


( 11 ) 


( 8 ) 

( 10 ) 
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The parameter r\/y can easily be eliminated, yielding. 


, _ . /I 2* 

(z— !=>) = r - V x +y 

or 

z = £> ± (R -\/c 2 +y 2 ) 

Another simple example is the boundary defined by: 

•n = £ 

where 

V = 1 
(9) reads: 

(x-£) +(y--n) = 0 

Eliminating | : 

x +y - 2r| = 0 
or 

_x+y 


(8) becomes: 


(z-t) 2 =(x-^) 2 + ,y -^, 2 

1. ,2 

= 2<x-y) 


or 


z = ^ ± — (x-y) 


When the boundary is defined by: 



i. e. : 


V = 0 

one immediately obtains from (9): 
x-i = o 

i. e. : 

(z-4) 2 = (y-b) 2 


z = i ± (y-b) 

Similarly for a boundary defined as: 
£ = a 

the solution is: 

z = ?, ± (x-a) 

The transformation of z into w is: 



.kw . , -1 
( — ) - smh 
c 



(ID 


( 9 ) 


( 10 ) 


and, for w Q = 0, 

w = sinh[k(z-C)J 

For the case of a rectangular plate. Dr. Naumann' s results presented 
above thus agree completely with the ones found in Chapter (B. 2): for instance, 
in region (4) of the rectangular domain (figure 3), for the boundary t] = b, we 
get: 

z-t, = ± (y-b) 
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i. e. : 


w = j sinh[±k(y-b)] = ^ sinh[k(y-b)], e = ± 1 

which is identical with the result we found. 

Another interesting remark concerns the subject of which boundary controls 
what part of the domain. The z-functions are planes which intersect the xy plane 
at an angle of 45°. Thus, the intersecting boundaries within the domain are simply 
the bissectors of the angles of the boundaries, which are the dotted lines 
represented in figure 3. The solution z can be visualized as being a roof over the 
boundary base with a pitch angle of 45° all the way around (as well as its exact 
complement, i. e. , a pitch angle of -45°). 

One can also visualize the function z to be like a sand hill on top of any 
shaped boundary with the sand of such a consistency that it can only support a 
pitch angle of 45°. 

To generalize the above results, it is obvious that the function z-£ is 
± p, if p is the shortest distance of the particular point in the xy plane to the 
given boundary 9D, which leads to the general solution for (1) given by: 

w = “|r sinh(kp) 

if the boundary condition is w = 0 along 9D. 

In case that w varies along 9D, i. e. , d£> £ 0, equation (9) is still 
applicable and reads: 

(z— C) jjj| = (x-£) + (y-r) )r| ' (9a) 

Again, equations (8) , (9a) , (10) , (11) , together with: 

jj|=f(£.*l) (12) 

uniquely define the function z=z(x,y) solution of (1), and the above procedure 
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can thus be used even if w varies in a known manner along the boundary 3D. 
Our alternate approach of Chapter (B. 2) could also have been extended to this 
more general ease, which is however of a limited practical value. 
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a o o n o o o o o o o o o n o r> r> r> r> 


APPENDIX 3 


A COMPUTER PROGRAM TO FIND THE OPTIMAL THICKNESS 
DISTRIBUTION OF A SANDWICH PLATE 


* FUNCTION M * 


CLASSICAL OPERATORS FOR THE PARTIAL DERIVATIVES 
STARTING FUNCTION IS W = S IN( PI*X/10 ) *S IN ( PI*Y/10 ) 

TEST PERFORMED ON THE DOMAIN l 0, 10} *{ 0, 10) 

STEP .2 ( 49*49 MESH POINTS) 

VALUE AT EACH POINT GOT BY SOLVING AW2+BW+C=0 
ONLY POSITIVE ROOT KEPT AS RESULT , STOP OTHERWISE 

REAL *8 W( 51 ,51 ) ,T( 51, 51) > 

1WXX,WXY >WYY> 

2WXXX > WXXY ,WXYY»WYYY, 

3WXXXX tWXXYY » WYYYY > 

4W2XXI 51,51) ,w2XY( 51,51 ) ,W2YYI 51, 51) , 

5W3XXX ( 51 > 51 ) * W3XXY ( 51 > 51 ) « W3XYY{ 5 It 51) , W3YYYC 5 1 , 5 1) > 
6W4XXXX(51,51), W4XXYY { 5 1 >51) ,W4YYYY< 51, 51) , 
7H»H2,H3,H4,A,B,C,DELTA,DIF2 > BETA2»HX» HY,TP,VIT, 

8AB1,AB2 

SETTING UP A FEW CONSTANTS 

BETA2=0.0239 
H=0«2 
H2=H*H 
H3=H2*H 
H4=H3*H 
N=49 
Nl=51 

A=-BETA2+10.64/H4 
K=1 

INITIALISATION 

DO 3 1=1, M 
DO 3 J=1 , N1 

W(I ,J )=DSIN( ( 1-1 )*H*3. 1416/ 10. ) *DSIN{ ( J-i)*H*3. 1416/10. ) 

3 CONTINUE 
C 

50 DO 1 I Y=1 »N 
DO 1 I X=1 , N 

B=-5.32/H4*lW{ I Y+l , IX+2J+WI IY+i,IX)+W( I Y+2 , IX+ 1 ) *W( IY,1X+1)) 

C 

C=1 . 34/16. 0/H4*( W( IY+2, IX+2 J + W ( I Y , I X)- W ( I Y+2 , 1 X) -W ( I Y, I X+2 )) ** 2 
1*(W( lY-tl, IX)+W( IY+1, IX+2) ) *#2/H4 
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no ooooonooo 


2+ ( W { I Y+2 »IX+1)+W( IY,IX+1) )**2/H4 

3+0.66*(W(IY+l,IX)?-W(IY+l,IX+2) ) *( W( I Y+2 , 1 X+1)+WU Y, IX+1 ) ) /H4-1 .0 
DELTA=B*B-4.0*A*C 
IF (OELTA.LT. 0.0) GOTO 30 
C 

10 W ( IY+1 , IX+1 )= (-B+DSQRT (DELTA) ) / ( 2.0*A) 

1 CONTINUE 
C 

DIF2=0.0 

HR! TE ( 6 »199 ) K 

K=K + 1 

IF (K.EQ.51 ) GOTO 30 
GOTO 50 

30 WRITE (6, 100) 

CALL FLEVELIK, -100. 0,2. 0,15, 0.015) 

00 4 J=l,26 
AB1=( J-1)*H 
WRITE (6,109) AB1 
DO 4 1=1,26 
AB2 = ( I— 1 ) *H 

WRITE(6,1C8) A82,W(J,I) 

4 CONTINUE 


* #4 ***** ***** Jit*####*## 

**** FUNCTION T ********************* 

* * * 35c * # * 5* * * * * * * * * * * * * * * * * ❖ # $ * * * $ A # # * 

INITIALISATION OF T 


DO 5 J~li N1 
DO 5 1=1, N1 

T ( I ,J )=DABS(1.5*DSIN( ( J-l ) *H*3 . 1416/5. ) *DSIN( ( 1-1) * H*3 . 1416/ 5. ) ) 
5 CONTINUE 


DO 77 J=1,N 
DO 77 1=1, N 

WXX=( W(J+1, 1+2) +W( J+l , I )— 2.0*W( J+l, 1+1) )/H2 
KYY=( h( J+2, I+l)+W< J,I+1)-2.0#W( J+l, 1+1) )/H2 
WXY— ( W( J+2 , I +2 ) +W ( J , I ) -W( J+2, I )-W ( J , I +2 ) ) /4. 0/H2 
WXXY= ( W( J+2 ,1 ) +W ( J+2, I +2 ) -W ( J , I )-W( J, I.+2 ) 

1—2. 0*W( J+2, I+1)+2.0*W(J»1+1) )/2.0/H3 
KXYY= ( W( J+2 , I +2 )— W ( J+2 , I ) +W ( J, I +2 ) — W( J , I ) 

1-2.0*W( J + l, I+2)+2-0*W< J+l ,1 ) )/2.0/H3 
h'XXYY= ( W ( J + 2 , 1 ) +W( J +2 , I +2 ) +W ( J , I ) +W ( J, I +2 ) 

1 - 2 . 0 * (W( J+l , I)+W(J+1» 1+2) +W(J, I.+1)+W(J+2,I+1))+4.0*W( J + l, I + 1D/H4 
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non o r> n 


IF { ( J .EQ.l) .AND. (I .EQ. 1)J GOTO 60 
IF ((J .EQ. N) .AND. (I .EQ. 1)1 GOTO 62 
IF ((J.EQ.N).AND.(I.EQ.N)) GOTO 64 
1F(( J.EQ.l) .AND.(I.EQ.N)) GOTO 66 
IF (I .EQ.l) GOTO 63 
IF (I-EQ.N) GOTO 70 
IF ( J .EQ. 1) GOTO 72 
IF (J.EQ.N) GOTO 74 

CASE OF AN INTERIOR POINT 

WXXXX=(U( J+l, I + 3)-4.0«[Jtl, I«-2)+W( J+ljI-1) 

1-4. 0*W { J+l , I ) +6.0*W { J+l , I+l ) )/H4 
WYYYY=(W( J+3, I+1)-4.0*W(J+2,1+1)+W(J-1, I + l) 
l-4.0*’ rt (J, H-l)+6.0*H(Jtl,Ul) )/H4 
WXXX=(W( J + l, I+3)-W< J+l, I-l)-2.0 : «W(J + l,I+2) + 2.()*W( J + l, I ) J/2.0/H3 
WYYY= ( W( J+3 » I + 1 )— W{ J-l , I+l)— 2»0*M( J+2 .1+1 )+2.0*W( J, I+l) I/2.0/H3 
GOTO 76 

CASE OF AN EDGE POINT 

60 WXXXX=(W( J+l , I+3J+5.0*W<J+1,1+1)-4.0*W( J+l, I+2)-2.0*W( J+l, I) )/H4 
1+(W(J+1, I+3)-3.0*W(J+l, 1+1) +2« Q*N( J+l , I) 3/3.0/H4 
WYYYY- { W( J + 3 1 1 + 1 ) +5 .O^W (J+l, I+l >-4. O^W (J+2,I+l)-2»0*W(JtI+l) >/H4 
1+ ( W( J +3 , I +1 )— 3 . Q*W ( J+l , I+l ) +2 .0*W( J » I + l ) ) /3.0/H4 
WXXX=(Vi( J + l, I*3)+W( J+l, I+l)-2.0*W(J+l,I+2) J/2.0/H3 
1-(W< J+l,I+3)-3.0*W(J+l, I+1)+2.0*W( J+l, I) )/6.0/H3 
WYYY=(W( J+3, I+l)+W( J+l, I+l)~2.0*W(J+2,I+l))/2.0/H3 
1— ( W ( J +3 , 1 +1 ) — 3 . 0*W ( J+l , I+l) +2.0*W( J» I+l) J/6.0/H3 
GOTO 76 

62 WXXXX=(W( J+l, I+3)+5.0*W(J+l,I+l)-4.0*W( J+l, 1+2 ) -2 . 0*W< J+l , I ) )/H4 
1+(W( J+l, I+3)-3.0*W(J+l,I+l)+2.0*W( J+l, I ) )/3.0/H4 
WYYYY = (W( J-l, I + l)+5.0<'W(J + l,I + l)i4.0*W( J, I+l ) -2 . ( J +2 , I+l ) 1/H4 
1+(W( J-l»I+l)-3.0*H{J+l,I + l)+2.0*VUJ+2,I+in/3.0/H4 
WXXX=(W( J + l, I+3)+W( J+1,I+1)-2.0*W(J+1, 1.+2))/2.0/H3 
1— < W C J+l,I+3)-3.0*W( J+l, I+1)+2.0*W( J+l, I ))/6.0/H3 
WYYY— — (W( J-l, I+l) +W (J+l, I+l ) -2.0*W( J , I +1 ) ) /2.0/H3 
l+(W(J-l,I+l)-3. 0*W( J+l, I+l)+2.0*W(J+2, I+l) I/6.0/H3 
GOTO 76 

64 WXXXX = (W< J+1,I-1)+5.0*W(J+1,I+1)-4.0*W( J+l , L) -2 .0*W { J+l , I +2 ) )/H4 

l+(W(J+lt I-1)-3.0*WC J+1,I+1)+2.0*W( J+l, 1+2) I/3.0/H4 
HYYYY=(W( J-l,I+l)+5.0*W(J+l,I+l)-4c0*W( J , I +1 ) -2 . 0*W ( J+2 , I+l ) )/H4 
1+(W( J-1,I+1)-3.0*H( J+1,I+1)+2.0*W(J+2,I+1 ) J/3.0/H4 
WXXX=-(W( J+l , I— 1 ) +W ( J+l , I+l )— 2.0*W( J+l , I) )/2.0/H3 
1+(H( J+l, I— 1) — 3.0 -W (J + l, I+l) +2 .0*W( J+l, I.+2) J/6.0/H3 
KYYY=-(W( J-l, I + l) +vf(J + l,I+l)-2.0*W( J,I+1) )/2.0/H3 
1+(W( J-l, I+1)-3.0*W( J+l, I+1)+Z.0*W(J+2,I+1) J/6.0/H3 
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GOTO 76 


C 

66 WXXXX=(W( J+l, 1-11+5. 0*W ( J+l, 1+1 > -4. 0*W{ J+l, I ) -2 .0*W { J +1 , 1+2) )/H4 

1+<W( J.+ l , I-1)-3.0*W(J+1, I+1)+2.0*K( J+l, 1+2) )/3.0/H4 
WYYYY=(W( J+3, I+1)+5.0*W( J+l,I+l)-4. 0*vl( J + 2, I + 1 ) -2 . 0*W< J , 1 + 1 ) )/H4 
l+(W(J+3, I+i)-3.0*W(J+l, I + 1)+2.0*'J(J,I + 1) ) /3.0/H4 
WXXX^-(W( JH , I-l)+W( J+1»I+1)-2.0*W< J + l , I ) )/2.0/H3 
1+(W( J+1,1-1 )-3.0*W( J+l, I+1)+2.0*W{ J+l, 1+2) J/6.0/H3 
WYYY={K{J+3,I+1)+W{ J+l, 1 + 1) -2. 0*WU +2,1+1) )/2.0/H3 
1-(W( J+3, I+1)-3.0*W( J+1,I+1)+2.0*W[ J,I+l))/6.0/H3 
GOTO 76 
C 

68 WXXXX=(W( J + l, I+3)+5.0*W( J + l, 1 + 1 )-.4.0*W( J + l, I+2)-2.0*W( J+l, I ) )/H4 

1+(W( J+l,I+3)-3.0*W{J+l,l+l)+2.0*W{J+l,I))/3.0/H4 
WYYYY= (W(J + 3»I+1 )— 4 • 0 J ! C W(J + 2»I+1)+W( J-1,1 + 1) 

1-4.0*W(J, I+1)+6.0*W(J+1, 1+1) )/H4 
WXXX=(W( J+l ,I+3)+W( J+l# I+1)-2.0*W(J+1, 1+2) )/2.0/H3 
1-(W( J+l,I+3)-3.0*W< J + l, I+1)+2.0*W( J + l, I) )/fe.0/H3 
WYYY=(W{ J+3,I+1)-W(J-1» I+l)-2.0*W(J+2, I+1)+2.0*W( J, 1 + 1) )/2.0/H3 
GOTO 76 
C 

7C WXXXX=(W( J+l,I-l)+5.0*W<J+l,I+l)-4. 0*K( J+l,I)-2.0*W(J+l,I+2) >/H4 
1+(W( J+l,i-l)-3.0*W< J+l, 1+1) +2.0*W(J+1, 1+2) ) /3 . 0/H4 
WYYYY=(W( J+3,I+1 )-4.0*W<J +2,1+1 )+W{ J-1,1+1) 

1-4* 0*W{ J, 1+1 ) +6.0*W( J+l , I +1) ) /H4 
WXXX=-(W{ J+l, I-l)+K(J+l,I+l)-2.0*W( J+l, I.) ) /2.0/H3 
1+{W( J+1,1-1 )-3.0*W< J+l, 1+1) +2. 0*W{ J+l, 1+2) )/6.0/H3 
WYYY=(W( J+3,I+1)-W( J-l,I+l)-2.0*W(J+2, I+1)+2.0*W( J, 1 + 1) J/2.0/H3 
GOTO 76 
C 

72 WXXXX=(W( J+l, I+3)-4.0*W(J+l,I>2)+W( J+l, I-i) 

1— 4.0 :< ‘V«(J + 1,I)+6.0*W(J + 1»I+1) )/H4 
WYYYY= { J+3 , 1+1 )+5 .0*W (J+1,I+1) —4. 0< C W(J+2»J + 1)-2.0 : <‘W(J,I+1) 1/H4 
1+(W( J+3,I+1)-3.0*WU + 1, I+1)+2.0*W{J, 1+1) J/3.0/H4 
WXXX=(W( J+l,I+3)-W(J+l» I-l)-2.0*Wl J+l,I+2)+2.0*W( J + l, I) 1/2.0/H3 
WYYY=('W(J+3, I+1J+W(J+1, 1+ 1) — 2. 0*W (J+2,11+1) 1/2.0/H3 
1— < W £ J+3,I+1 )-3.0*W< J+l, I+1)+2.0*W( J , 1+1) )/6.0/H3 
GOTO 76 
C 

74 WXXXX=(W( J+l, I+3)-4.0*W( J+l,I+2)+W( J+l, 1-1) 
1-4.0*WCJ+1,I)+6.0*W(J+1,I+1))/H4 
WYYYY = ( W( J-L, I+l)+5.0*W<J+l,I+l)-4.0*W{ J, I +1 > -2. ( J+2 , I +1 ) )/H4 
1 + (W< J-1,1+1 )-3.0*W< J+l, I+1)+2.0*W(J+2,I+1) J/3.0/H4 
WXXX=(W( Jtl,I+3)-W( J+l, I-1)-2.0*W(J+1, I+2)+2.0*WC J+l, I ) 1/2.0/H3 
WYYY— — ( n( J-l, I+1)+VHJ+1,I+1)-2.0*W{ J»I+l) )/2.0/H3 
1 + (W{J— 1,1+1 )— 3.0#V) (J+l, I + 1)+2.0*W< J+2,1 + 1) )/6.0/H3 

76 W2XX( J+l , 1+1 )=WXX 
W2XY( J+l, 1+ 1) =WXY 
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W2YY ( J+l, 1+1 )=WYY 
V13XXX { J+l, l+l)=WXXX 
W3XXY{ J+l ,I+1)=WXXY 
W3XYY (J+1,1+1 )=WXYY 
W3YYY ( J+l f I +1 } =WYYY 
W4XXXXI J+l, 1+1)=WXXXX 
W4XXYY< J+l, l+l)=WXXYY 
W4YYYY( J+l, I +1 )=WYYYY 
77 CONTINUE 
C 

K=1 

55 DIF2=0.0 

C 

DO 6 J = 1,N 
DO 6 1=1, N 
C 

A=-2.66/H2* (W2XXIJ+1, 1+1) +W2YYIJ +1,1+1) ) 

1+W4XXXX( J + l, I+1)+2.0*W4XXYY( J + l, 1+1) +W4YYYY( J + 1,1 + 1 ) 

2-BETA2*W( J+1,1+1) 

B=(W2XX( J+l, I+1)+0.33*W2YY(J+1, 1 + 1) ) *( T ( J+l, I.+2 ) +T( J+l , I ) J /H2 
1+1. 34# (T( J+2 , I+2I+TI J , I )— T( J+2, I )—T(J,I+2)) *W2XY (J+l,I+l)/4.0/H2 
2+(0.-33*W2XX( J +1,1+1 )+W2YY( J+1,1+1 ))*(T( J+2, 1+l ) +-T ( J , 1+1 ) } /H2 
2+(W3XXX( J+l ,1+1 )+W3XYY( J+1,1+1) )*(T< J+l, 1+2) -T( J+l, I))/H 
4+ C W3XXYI J +1 , 1+1 ) +W3YYY (J + 1,1+1) )*(T ( J+2,I+1)-T (J, I + 1 1 )/H 
5-0. 63 2#BETA2*W( J+1,1+1) 

C 

TP=-B/A 

DIF2=DIF2+(T(I+1,J+1J-TP)<‘*2 
T (J+1,1+1 J=TP 
6 CONTINUE 

C 

WRITE (6,200 ) K.DIF2 

IF( (K.EQ.150I .OR. ( DI F2. LT .0.001) ) GOTD 80 
K=K+1 
GOTO 55 
C 

C FINAL RESULTS FOR FUNCTION T 

C 


80 

WRITE 

(6,107) 



DO 7 

J=1 , 26 



AB1= ( 

J— 1 ) #H 



WRITE 

(6,110) 

AB1 


DO 7 

1=1,26 



AB2=( 

1-1) #H 



WRITE 

(6,111) 

AB2,T ( J, I ) 

7 

CONTI 

NUE 



c 

HX=10./49. 
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HY=HX 
VIT=0.0 
OD 2 1-2 , 49 
DO 2 J-=2 » 49 

VIT=VIT*T{J-1»I-1)+T( J+2t I-l)+T( J-1 jI+2>+ 

1 T ( J+2f J t I-ll+Tt J+li I-D + 

2 1 ( J-l»IH-T{J+2»I)f-T( J-l,Dl)tTlJ+2iIHI 

3 +T« J f I+2)+Tr J + l,H-2) H-169.0*(T( J,I )+T{ J+1»I ) 

4 +T{ J » 1*1 ) -*-T( J-*-l» I>1) ) 

2 CONTINUE 

C 

VI T= V l T*HX*HY/24. 0/24.0 
C 

WRITE ( 6 y 600 } VIT 
C 

CALL ELEVEL(T,-100.0,2.0,15,0.015) 

100 FORMAT (lHlt IOX» 'FUNCTION W ',///> 

107 FORMAT(1H1.10X, 'FUNCTION T *»//) 

1C8 FORMAT ( 1H , 15X, E16 .7 , 15X, El 6 .7 , / ) 

109 FORMAT ClHli • LINE Y=« , E16. 7 ,// ,25X, ' X • , 15X, » W( P) • , / ) 

110 FORMAT ( 1H1» ' L I NE Y=* ,E 16 .7, //, 25X, 'X* , 15X, » T ( P ) ’ » / ) 

111 FORMAT-! 1H , 15X, E16. 7 , 15X, E16. 7 , / ) 

199 FORMAT ( 1H ,' ITERATION N. :*»I3*/I 

200 FORMAT ( 1H , 'ITERATION N. ; * , I 3 , 5X , ' DELTA2 S',E16.7,/) 

600 FORMAT (IHlt 'VALUE OF THE INTEGRAL IS :',E16.7,//) 

RETURN 

END 
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SUBROUTINE FLEVEL(Z»WL,UL,NBL,TOL) 
REAL *8 Z »W L» UL , TOL 
DIMENSION Z ( 5 1, 51) 


C 

C TRACES TEE LEVEL-LINES FOR A TABULATED FUNCTION { VALUES IN Z(J,D) 
C WL IS THE LEVEL UNDER WHICH •?• IS PRINTED 
C UL SAME FUNCTION FOR THE UPPER LEVEL 
C N3L IS THE NUMBER OF LEVELS (STEP .1 ASSUMED) 

C TOL IS THE HALF OF THE WIDTH OF THE BOUND 
C 

C PROGRAM PRINTS OUT THE LEVELS IDENTIFIED BY LETTERS ON ONE PAGE 
C 

INTEGER LIT *2(15}/*A *,*B *,»C *,'0 *,*E *,*F ' , *G *,*H ', 

1 *1 J S'K • » ' L * » * M ' » 1 N *,'0 '/ 

INTEGER PAGE *2(51,51) 

INTEGER *2 QUEST/* ? «/, BLANK./* ’/ 

C 

DC 1 1=1,51 
DO 1 J=1 , 51 
PAGE ( J ,1 ) =BLANK 

1 CONTINUE 
C 

DO 2 J=2 , 50 
DO 2 1=2,50 

IF ( ( Z(J, I) .LT.WD.OR. ( Z( J , I ) .GT.UL) ) PAGE( J, I )=QUEST 

2 CONTINUE 
C 

DO 3 K=l, NBL 
VAL=0 » 1*K 
C 

DO 4 J=1 » 51 
DO 4 1=1,51 

IF ( ( Z(J, I) .GT. (VAL-TOL) ) » AND, (Z(J,I).LT. (VAL+TOL) ) ) 

1PAGEI J,I)=LIT(K) 

4 CONTINUE 

3 CONTINUE 
C 

DO 5 J=1 , 51 

WRITE (6,500) { PAGEl 52— J , I ) ,1=1,51) 

5 CGNTINUE 
C 

5C0 FORMAT (1H ,20X,51A2) 

WRI TE (6, 501 ) 

501 FORMAT(lHl) 

RETURN 

END 


149 



REFERENCES 


1. Shell, C. Y. and Prager, W. , "Recent Developments in Optimal Structural 
Design, " Applied Mechanics Review , Vol. 21, No. 10, October 1968, 

pp. 985-992. 

2. Haug, E. J. , Jr. , Streeter, T. D. and Newell, R. S. , "Optimal Design 
of Elastic Structural Elements, " SY-RI-69, April 1969, Systems Analysis 
Directorate, U. S. Army Weapons Command, Rock Island, Illinois. 

3. MacCart, B. R. , Haug, E. J. , Jr. and Streeter, T. D. , "Optimal Design 
of Structures with Constraints on Natural Frequency", AIAA Journal , 

Vol. 8, No. 6, June 1970, pp. 1012-1019. 

4. Fox, R. L. and Kapoor, M. P. , "Structural Optimization in the Dynamic 
Response Regime: A Computational Approach", AIAA Journal , Vol. 8, No. 

10, October 1970, pp. 1798-1804. 

5. Rubin, C. P. , "Minimum- Weight Design of Complex Structures Subject to a 
Frequency Constraint", AIAA Journal , Vol. 8, No. 5, May 1970, pp. 923-927. 

6. Turner, M. J. , "Optimization of Structures to Satisfy Flutter Requirements", 
AIAA Structural Dynamics and Aeroelasticity Specialist Conference, New 
Orleans, La., Apr. 16-17, 1969, Volume of Technical Papers , pp. 1-8. 

7. Turner, M. J. , "Design of Minimum-Mass Structures with Specified Natural 
Frequencies", AIAA Journal , Vol. 5, No. 3, March 1967, pp. 406-412. 

8. Keller, J. B. and Niordson, F. I. , "The Tallest Column", Journal of 
Mathematics and Mechanics , Vol. 16, No. 5, Nov. 1966, pp. 433-446. 

9. Niordson, F. I. , "On the Optimal Design of a Vibration Beam", Quarterly 
of Applied Mathematics , Vol. XXIII, No. 1, April 1965, pp. 47-53. 

10. Olhoff, N. , "Optimal Design of Vibrating Circular Plates", Department of 

Applied Mechanics, The Technical University of Denmark, Copenhagen, 1969. 


150 


11. Prager, W. and Taylor, J. E. , "Problems of Optimal Structural Design", 
Journal of Applied Mechanics, Transactions ASME , Vol. 35E, No. 1, 

March 1968, pp. 102-106. 

12. Taylor, J. E. , "Optimum Design of a Vibrating Bar with Specified Minimum 
Cross-Section", AIAA Journal , Vol. 6, No. 7, July 1968, pp. 1379-1381. 

13. Ashley, H. and McIntosh, S. C. , Jr. , "Application of Aeroelastic Constraints 
in Structural Optimization", Proceedings of the Twelfth International 
Congress of Applied Mechanics , Stanford University, Stanford, California, 
Aug. 1968, Springer-Verlag, New York. 

14. Armand, J. -L. and Vitte, W. J. , "Foundations of Aeroelastic Optimization 
and Some Applications to Continuous Systems", SUDAAJR 390, Jan. 1970, 
Department of Aeronautics and Astronautics , Stanford University, Stanford, 
California. 

15. Weisshaar, T. A. , "An Application of Control Theory Methods to the 
Optimization of Structures Having Dynamic or Aeroelastic Constraints", 
SUDAAJR 412, Oct. 1970, Department of Aeronautics and Astronautics, 
Stanford University, Stanford, California. 

16. McIntosh, S. C. , Jr. , Weisshaar, T. A. and Ashley, H. , "Progress in 
Aeroelastic Optimization- Analytical Versus Numerical Approaches", 

SUDAAR 383, July 1969, Department of Aeronautics and Astronautics, 
Stanford University, Stanford, California. 

17. Sippel, D. L. , "Minimum-Mass Design of Structural Elements and Multi- 
Element Systems with Specified Natural Frequencies", Ph. D. Dissertation, 
September 1970, University of Minnesota. 

18. Johnson, M. R. , "Optimum Frequency Design of Structural Elements", 

Ph. D. Dissertation, 1968, Department of Engineering Mechanics , The 
University of Iowa. 


151 



19. H aug, E. J.j Pan, K. C. ; and Streeter, T. D. , "A Computational Method 
for Optimal Structural Design I", Int. J. Num. Meth. Engr . , 1972* 

20. Butkovskiy, A. G. and Leraer, A. Ya. , "The Optimal Control of Systems 
with Distributed Parameters", Automation and Remote Control , Vol. 21, 

No. 6, Dec. 1960, pp. 472-477. 

21. Robinson, A. C. , "A Survey of Optimal Control of Distributed-Parameter 
Systems", ARL 69-0177, Nov. 1969, Office of Aerospace Research, U. S. 
Air Force. 

22. Butkovskiy, A. G. , Distributed Control Systems , American Elsevier 
Publishing Co. , Inc. , New York, 1969. 

23. Lions, J. L. , ContrSle Optimal de Systemes Gouvern^s par des Equations 
aux Derivees Partielles , Dunod et Gauthier-Villars, Editeurs, Paris, 1968. 

24. Lur'e, K. A. , "The Mayer-Bolza Problem for Multiple Integrals and the 
Optimization of the Performance of Systems with Distributed Parameters", 
Applied Mathematics and Mechanics (PMM) , Vol. 27, No. 5, March 1963, 
pp. 1284-1299. 

25. Lur’ e, K. A. , "The Mayer-Bolza Problem for Multiple Integrals: Some 
Optimum Problems for Elliptic Differential Equations Arising in Magneto- 
hydrodynamics", in Topics in Optimization , Chapter 5, pp. 147-198, 

G. Leitmann, Editor, Academic Press, New York, 1967. 

26. Lur'e, K. A. , "Optimum Control of Conductivity of a Fluid Moving in a 
Channel in a Magnetic Field", Journal of Applied Mathematics and 
Mechanics , Vol. 28, No. 2, 1964, pp. 316-327. 

27. Bryson, A. E. , Jr. and Ho, Y. -C. , Applied Optimal Control, The Blaisdell 
Publishing Co. , Waltham, Massachusetts, 1969. 

28. Pontryagin, L. S. , et al. , The Mathematical Theory of Optimal Processes , 
Interscience, New York, 1962. 


152 


29. Butkovskiy, A. G. , "The Maximum Principle for Optimum Systems with 
Distributed Parameters", Automation and Remote Control , Vol. 22, No. 10, 
March 1962, pp. 1156-1169. 

30. Egorov, A. I. , "Optimality Conditions of Systems Containing Components 
with Distributed Parameters", in Mathematical Theory of Control , A. V. 
Balakrishnan and Lucien W. Neustadt, Editors, Academic Press, New York, 
1967, pp. 322-329. 

31. Neustadt, L. W. , "A General Theory of Extremals", Journal of Computer 
and System Sciences , Vol. 3, Feb. 1969, pp. 57-92. 

32. Meirovitch, L. , Analytical Methods in Vibrations , The MacMillan Company, 
New York, 1967. 

33. Courant, R. and Hilbert, D. , Methods of Mathematical Physics , Vol. 1, 
Interscience Publishers, Inc. , New York, 1953. 

34. MacLachlan, N. W. , Bessel Functions for Engineers , 2nd ed. , Oxford 
University Press, New York, 1961. 

35. Garabedian, P. R. , Partial Differential Equations , John Wiley and Sons, 

Inc. , New York, 1964. 

36. Reissner, E. , "Small Bending and Stretching of Sandwich-Type Shells", 

NACA Report 975 . 

37. Reissner, E. , "OnBending of Elastic Plates", Quarterly of Applied 
Mathematics , Vol. V, No. 1, 1947, pp. 55-68. 

38. Hoff, N. J. , "Bending and Buckling of Rectangular Sandwich Plates", 

NACA Technical Note 2225 . 

39. Ames, W. F. , Nonlinear Partial Differential Equations in Engineering , 
Academic Press, New York, N. Y. , 1965. 

40. Forsyth, A. R. , Theory of Differential Equations , 4 parts in 6 volumes, 
Cambridge University Press, 1890-1906; also reprinted by Dover Publications, 
New York, N. Y. 1955. 


153 



41 . 


Ampere, A. M. , "Considerations Generates sur les Integrates des 
Equations aux DifferentieU.es Partielles", Journal de 1' Ecole Poly technique , 
XVIIIe Caliier, Tome X, Janvier 1815, pp. 549-611. 

42. Ampere, A. M. , "Memoire Contenant 1' Application de la Theorie Exposee 
dans le XVUe Cahier du Journal de 1' Ecole Polytechnique, a 1' Integration 
des Equations aux Differentielles Partielles du Premier et du Second 
Ordre", Journal de P Ecole Royale Polytechnique , XVIIIe Cahier, Tome XI, 
Janvier 1820, pp. 1-188. 

43. Timoshenko, S. P. and Woinowsky-Krieger , S. , Theory of Plates and 
Shells , 2nd edition, McGraw-Hill Book Company, New York, N. Y. , 1959. 

44. Reissner, E. , "Remark on the Theory of Bending of Plates of Variable 
Thickness", Journal of Mathematics and Physics , Vol. 16, 1937, pp. 43-45. 

45. Suhubi, E. , "Design of Plates for Minimum-Weight", Istanbul Teknik 
Universitesi Bulteni , Vol. 14, No. 1, 1961, pp. 11-30. 

46. Prager, W. , "Optimality Criteria Derived from Classical Extremum 
Principles", in An Introduction to Structural Optimization , M. Z. Cohn, 

Ed. , University of Waterloo, Waterloo, Ontario, Canada, 1969, pp. 165-178. 

47. Prager, W. , "Optimality Criteria in Structural Design", Proceedings of the 
National Academy of Sciences, Vol. 61, 1968, pp. 794-796. 

48. Munk, M. M. , "Isoperimetrische Probleme aus der Theorie des Fluges", 
Gottingen Dissertation, 1918. 

49. Forsythe, G. E. and Wasow, W. R. , Finite-Difference Methods for Partial 
Differential Equations , John Wiley and Sons, Inc. , New York, N. Y. 1960. 

50. Chattopadhyay, R. , "Control of Distributed-Parameter Systems", 66-32, 

May 1966, Department of Engineering, University of California, Los Angeles, 
California. 


154 


CR-2044 


— 32 


